zhyang-dev/LIGGGHTS-PFM
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

particle fragmentation: initial commit

Browse Source
This commit is contained in:
Daniel
2014-11-19 00:05:08 +01:00
parent ea1e50a33f
commit b4808232a9
14 changed files with 2417 additions and 0 deletions
Download Patch File Download Diff File Expand all files Collapse all files

824
src/fix_break_particle.cpp Normal file
View File

@ -0,0 +1,824 @@
/* ----------------------------------------------------------------------
LIGGGHTS - LAMMPS Improved for General Granular and Granular Heat
Transfer Simulations
LIGGGHTS is part of the CFDEMproject
www.liggghts.com | www.cfdem.com
Christoph Kloss, christoph.kloss@cfdem.com
Copyright 2009-2012 JKU Linz
Copyright 2012- DCS Computing GmbH, Linz
LIGGGHTS is based on LAMMPS
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
This software is distributed under the GNU General Public License.
See the README file in the top-level directory.
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
Contributing authors:
Christoph Kloss (JKU Linz, DCS Computing GmbH, Linz)
Daniel Queteschiner (JKU Linz)
------------------------------------------------------------------------- */
#include "math.h"
#include "stdlib.h"
#include "string.h"
#include "fix_break_particle.h"
#include "atom.h"
#include "atom_vec.h"
#include "update.h"
#include "comm.h"
#include "modify.h"
#include "vector_liggghts.h"
#include "mpi_liggghts.h"
#include "domain.h"
#include "random_park.h"
#include "memory.h"
#include "error.h"
#include "fix_property_atom.h"
#include "fix_particledistribution_discrete.h"
#include "fix_template_fragments.h"
#include "particleToInsert_fragments.h"
#include "pair_gran.h"
#include "force.h"
#include "neigh_list.h"
#include <string>
#include <algorithm>
#include "math_const.h"
using namespace MathConst;
using namespace LAMMPS_NS;
using namespace FixConst;
enum{BC_ENERGY, BC_FORCE, BC_VON_MISES};
#define LMP_DEBUGMODE_FIX_BREAK_PARTICLE true
#define LMP_DEBUG_OUT_FIX_BREAK_PARTICLE screen
/* ---------------------------------------------------------------------- */
FixBreakParticle::FixBreakParticle(LAMMPS *lmp, int narg, char **arg) :
FixInsert(lmp, narg, arg)
{
// set defaults first, then parse args
init_defaults();
bool hasargs = true;
while (iarg < narg && hasargs) {
hasargs = false;
if (strcmp(arg[iarg],"fMat") == 0) {
if (iarg+2 > narg) error->fix_error(FLERR,this,"not enough arguments: missing fMat value");
fMat = atof(arg[iarg+1]);
if(fMat <= 0. ) error->fix_error(FLERR,this,"'fMat' must be > 0");
iarg += 2;
hasargs = true;
} else if (strcmp(arg[iarg],"energy_threshold") == 0) {
if (iarg+2 > narg) error->fix_error(FLERR,this,"not enough arguments: missing energy_threshold value");
threshold = atof(arg[iarg+1]);
if (threshold < 0.) error->fix_error(FLERR,this,"'energy_threshold' must be >= 0");
breakage_criterion = BC_ENERGY;
iarg += 2;
hasargs = true;
} else if (strcmp(arg[iarg],"force_threshold") == 0) {
if (iarg+2 > narg) error->fix_error(FLERR,this,"not enough arguments: missing force_threshold value");
threshold = atof(arg[iarg+1]);
if (threshold <= 0.) error->fix_error(FLERR,this,"'force_threshold' must be > 0");
breakage_criterion = BC_FORCE;
iarg += 2;
hasargs = true;
} else if (strcmp(arg[iarg],"von_mises_stress") == 0) {
if (iarg+2 > narg) error->fix_error(FLERR,this,"not enough arguments: missing von_mises_stress value");
threshold = atof(arg[iarg+1]);
if (threshold <= 0.) error->fix_error(FLERR,this,"'von_mises_stress' must be > 0");
breakage_criterion = BC_VON_MISES;
iarg += 2;
hasargs = true;
} else if (strcmp(arg[iarg],"min_radius") == 0) {
if (iarg+2 > narg) error->fix_error(FLERR,this,"not enough arguments: missing min_rad value");
min_break_rad = atof(arg[iarg+1]);
if (min_break_rad < 0.) error->fix_error(FLERR,this,"'min_break_rad' must be >= 0");
iarg += 2;
hasargs = true;
} else {
std::string unknown("unknown keyword ");
unknown += arg[iarg];
error->fix_error(FLERR,this,unknown.c_str());
}
}
// no fixed insertion target (such as # particles) since insertion triggered by break-ups
ninsert_exists = 0;
// turn off overlap check and turn off start stats since both not required
check_ol_flag = 0;
print_stats_start_flag = 0;
breakdata = NULL;
fix_break = NULL;
fix_breaker = NULL;
fix_collision_factor = NULL;
tag_max = 0;
// execute end of step
nevery = 1;
n_break = 0;
mass_break = 0.;
/////if(maxrad >= 1.) error->fix_error(FLERR,this,"Particle distribution must be relative, max radius mus be < 1");
}
/* ---------------------------------------------------------------------- */
void FixBreakParticle::post_create()
{
if (!fix_break) {
char * breakvar_name = new char[7+strlen(id)];
sprintf(breakvar_name,"break_%s",id);
const char * fixarg[9];
fixarg[0] = breakvar_name;
fixarg[1] = "all";
fixarg[2] = "property/atom";
fixarg[3] = breakvar_name;
fixarg[4] = "scalar"; //NP 1 scalar per particle to be registered
fixarg[5] = "yes"; //NP restart yes
fixarg[6] = "yes"; //NP communicate ghost yes
fixarg[7] = "no"; //NP communicate rev no
fixarg[8] = "0.";
fix_break = modify->add_fix_property_atom(9,const_cast<char**>(fixarg),style);
delete [] breakvar_name;
}
if (!fix_breaker) {
char * breakvar_name = new char[9+strlen(id)];
sprintf(breakvar_name,"breaker_%s",id);
const char * fixarg[9];
fixarg[0] = breakvar_name;
fixarg[1] = "all";
fixarg[2] = "property/atom";
fixarg[3] = breakvar_name;
fixarg[4] = "scalar"; //NP 1 scalar per particle to be registered
fixarg[5] = "yes"; //NP restart yes
fixarg[6] = "yes"; //NP communicate ghost yes
fixarg[7] = "no"; //NP communicate rev no
fixarg[8] = "0.";
fix_breaker = modify->add_fix_property_atom(9,const_cast<char**>(fixarg),style);
delete [] breakvar_name;
}
if (!fix_collision_factor) {
char * breakvar_name = new char[10+strlen(id)];
sprintf(breakvar_name,"break_cf_%s",id);
const char * fixarg[9];
fixarg[0] = breakvar_name;
fixarg[1] = "all";
fixarg[2] = "property/atom";
fixarg[3] = breakvar_name;
fixarg[4] = "scalar"; //NP 1 scalar per particle to be registered
fixarg[5] = "yes"; //NP restart yes
fixarg[6] = "yes"; //NP communicate ghost yes
fixarg[7] = "no"; //NP communicate rev no
fixarg[8] = "1.";
fix_collision_factor = modify->add_fix_property_atom(9,const_cast<char**>(fixarg),style);
delete [] breakvar_name;
}
}
/* ---------------------------------------------------------------------- */
void FixBreakParticle::pre_delete(bool unfixflag)
{
modify->delete_fix(fix_break->id);
}
/* ---------------------------------------------------------------------- */
FixBreakParticle::~FixBreakParticle()
{
if(breakdata) memory->destroy(breakdata);
}
/* ---------------------------------------------------------------------- */
void FixBreakParticle::init_defaults()
{
threshold = -1.0;
fMat = 1.0;
breakage_criterion = BC_ENERGY;
fix_fragments = NULL;
min_break_rad = 0.0;
}
/* ---------------------------------------------------------------------- */
void FixBreakParticle::init()
{
FixInsert::init();
pair_gran = static_cast<PairGran*>(force->pair_match("gran", 0));
dnum = pair_gran->dnum();
deltaMaxOffset = pair_gran->get_history_value_offset("deltaMax");
siblingOffset = pair_gran->get_history_value_offset("sibling");
collisionFactorOffset = pair_gran->get_history_value_offset("collisionFactor");
impactEnergyOffset = pair_gran->get_history_value_offset("impactEnergy");
if (deltaMaxOffset < 0 || impactEnergyOffset < 0)
error->fix_error(FLERR,this,"failed to find deltaMaxOffset/impactEnergy value in contact history");
}
/* ----------------------------------------------------------------------
calculate ninsert, insert_every, ninsert_per, massinsert, flowrates etc
also perform error checks
------------------------------------------------------------------------- */
void FixBreakParticle::calc_insertion_properties()
{
// error checks
if (threshold < 0.)
error->fix_error(FLERR,this,"you have to specify a threshold value");
if (nflowrate > 0. || massflowrate > 0.)
error->fix_error(FLERR,this,"specifying 'nflowrate' or 'massflowrate' is not allowed");
if (ninsert > 0 || massinsert > 0.)
error->fix_error(FLERR,this,"specifying 'nparticles' or 'mass' is not allowed");
if (insert_every <= 0)
error->fix_error(FLERR,this,"specifying 'every' must be > 0");
// fix holding particle fragments
if(fix_distribution->n_particletemplates() != 1)
error->fix_error(FLERR,this,"fix of type particledistribution/discrete must hold exactly one template");
if(strcmp(fix_distribution->particletemplates()[0]->style,"particletemplate/fragments"))
error->fix_error(FLERR,this,"fix of type particledistribution/discrete must hold exactly one template of type fix particletemplate/fragments");
fix_fragments = static_cast<FixTemplateFragments*>(fix_distribution->particletemplates()[0]);
// do not need ninsert_per
}
/* ---------------------------------------------------------------------- */
int FixBreakParticle::setmask()
{
int mask = FixInsert::setmask();
mask |= PRE_FORCE;
mask |= END_OF_STEP;
return mask;
}
/* ---------------------------------------------------------------------- */
double FixBreakParticle::insertion_fraction()
{
if(n_break_this == 0) return 0.0;
return static_cast<double>(n_break_this_local)/static_cast<double>(n_break_this);
}
/* ---------------------------------------------------------------------- */
inline int FixBreakParticle::is_nearby(int i)
{
// need not check overlap with existing particles since we
// know space originally taken by deleted particles is free
return 0;
}
/* ---------------------------------------------------------------------- */
void FixBreakParticle::pre_force(int)
{
// set sibling contact flags and collision factor in contact_history
if(n_break_this_local > 0) {
double **x = atom->x;
double *radius = atom->radius;
int *mask = atom->mask;
int inum = pair_gran->list->inum;
int * ilist = pair_gran->list->ilist;
int * numneigh = pair_gran->list->numneigh;
int ** firstneigh = pair_gran->list->firstneigh;
double ** firsthist = pair_gran->listgranhistory->firstdouble;
for (int ii = 0; ii < inum; ++ii) {
const int i = ilist[ii];
if (mask[i] & groupbit) {
const double xtmp = x[i][0];
const double ytmp = x[i][1];
const double ztmp = x[i][2];
const double radi = radius[i];
double * const allhist = firsthist[i];
int * const jlist = firstneigh[i];
const int jnum = numneigh[i];
for (int jj = 0; jj < jnum; ++jj) {
const int j = jlist[jj];
const double delx = xtmp - x[j][0];
const double dely = ytmp - x[j][1];
const double delz = ztmp - x[j][2];
const double rsq = delx * delx + dely * dely + delz * delz;
const double radj = radius[j];
const double radsum = radi + radj;
if (rsq < radsum * radsum) {
double * contact_history = &allhist[dnum*jj];
if (contact_history[deltaMaxOffset] == 0.0) {
contact_history[siblingOffset] = 1.0;
contact_history[collisionFactorOffset] = fix_collision_factor->vector_atom[i];
}
}
}
}
}
}
}
/* ---------------------------------------------------------------------- */
void FixBreakParticle::end_of_step()
{
double *flag = fix_break->vector_atom;
double *breaker = fix_breaker->vector_atom;
int nlocal = atom->nlocal;
int *mask = atom->mask;
int *tag = atom->tag;
double *radius = atom->radius;
int inum = pair_gran->list->inum;
int * ilist = pair_gran->list->ilist;
int * numneigh = pair_gran->list->numneigh;
int ** firstneigh = pair_gran->list->firstneigh;
double ** firsthist = pair_gran->listgranhistory->firstdouble;
switch (breakage_criterion) {
case BC_ENERGY:
for (int ii = 0; ii < inum; ++ii) {
const int i = ilist[ii];
if (mask[i] & groupbit && radius[i] > min_break_rad) {
double * const allhist = firsthist[i];
int * const jlist = firstneigh[i];
const int jnum = numneigh[i];
double breaker_energy = 0.0;
int breaker_tag = 0;
for (int jj = 0; jj < jnum; ++jj) {
const int j = jlist[jj];
double * contact_history = &allhist[dnum*jj];
const double impact_energy_limited = std::max(0.0, contact_history[impactEnergyOffset] - 0.5*threshold/radius[i]);
if(breaker_energy < impact_energy_limited) {
breaker_energy = impact_energy_limited;
breaker_tag = tag[j];
}
flag[i] += impact_energy_limited;
}
if(breaker_energy > 0.0) {
double probability = 1.0 - exp(-fMat * 2.0*radius[i] * flag[i]);
if (probability > random->uniform()) {
breaker[i] = static_cast<double>(breaker_tag);
}
}
//if(i < 4 && flag[i] != 0.0)
//fprintf(screen,"FixBreakParticle impact energy %d: %f\n", i, flag[i]);
}
}
break;
case BC_FORCE:
for(int i = 0; i < nlocal; ++i) {
if (mask[i] & groupbit && radius[i] > min_break_rad) {
/// TODO
}
}
break;
case BC_VON_MISES:
/// TODO
break;
default:
break;
}
}
/* ---------------------------------------------------------------------- */
void FixBreakParticle::pre_insert()
{
int i,ibreak;
int nlocal = atom->nlocal;
int *mask = atom->mask;
double **x = atom->x;
double **v = atom->v;
double *radius = atom->radius;
double *rmass = atom->rmass;
double *flag = fix_break->vector_atom;
double *breaker = fix_breaker->vector_atom;
AtomVec *avec = atom->avec;
// count # of particles to break
n_break_this_local = 0;
mass_break_this_local = 0.;
int inum = pair_gran->list->inum;
int * ilist = pair_gran->list->ilist;
int * numneigh = pair_gran->list->numneigh;
int ** firstneigh = pair_gran->list->firstneigh;
double ** firsthist = pair_gran->listgranhistory->firstdouble;
for (int ii = 0; ii < inum; ++ii) {
const int i = ilist[ii];
if (mask[i] & groupbit && breaker[i] > 0) {
double * const allhist = firsthist[i];
int * const jlist = firstneigh[i];
const int jnum = numneigh[i];
bool found_breaker = false;
for (int jj = 0; jj < jnum; ++jj) {
const int j = jlist[jj];
if (atom->tag[j] == breaker[i]) {
found_breaker = true;
double * contact_history = &allhist[dnum*jj];
if (contact_history[deltaMaxOffset] < 0.0) {
++n_break_this_local;
mass_break_this_local += rmass[i];
flag[i] = -(1.0 - exp(-fMat * 2.0*radius[i] * flag[i])); // sign indicates breakage
break;
}
}
}
if(!found_breaker) {
fprintf(screen,"FixBreakParticle::pre_insert: breaker not found!\n");
}
}
}
// tally stats
MPI_Sum_Scalar(n_break_this_local,n_break_this,world);
n_break += n_break_this;
MPI_Sum_Scalar(mass_break_this_local,mass_break_this,world);
mass_break += mass_break_this;
// virtual contact force calculation to restore overlap energy of contacting particles
delta_v.clear();
std::multimap<int,int> processed_pairs;
for (int ii = 0; ii < inum; ++ii) {
const int i = ilist[ii];
if (mask[i] & groupbit && flag[i] < 0.0) { // this one breaks
int * const jlist = firstneigh[i];
const int jnum = numneigh[i];
{
std::map<int, std::vector<double> >::iterator it = delta_v.find(i);
if (it == delta_v.end()) {
std::vector<double> delta_v0(3,0.0);
delta_v[i] = delta_v0;
}
}
for (int jj = 0; jj < jnum; ++jj) {
const int j = jlist[jj];
// check if this pair has already been processed
bool processed = false;
std::pair <std::multimap<int,int>::iterator, std::multimap<int,int>::iterator> ret;
ret = processed_pairs.equal_range(j);
for (std::multimap<int,int>::iterator it=ret.first; it!=ret.second; ++it) {
if(it->second == i) {
processed = true;
break;
}
}
if (processed) {
continue;
}
virtual_x_i.clear();
virtual_x_i.push_back(x[i][0]);
virtual_x_i.push_back(x[i][1]);
virtual_x_i.push_back(x[i][2]);
virtual_v_i.clear();
virtual_v_i.push_back(v[i][0]);
virtual_v_i.push_back(v[i][1]);
virtual_v_i.push_back(v[i][2]);
virtual_x_j.clear();
virtual_x_j.push_back(x[j][0]);
virtual_x_j.push_back(x[j][1]);
virtual_x_j.push_back(x[j][2]);
virtual_v_j.clear();
virtual_v_j.push_back(v[j][0]);
virtual_v_j.push_back(v[j][1]);
virtual_v_j.push_back(v[j][2]);
virtual_f_ij.clear();
virtual_f_ij.push_back(0.0);
virtual_f_ij.push_back(0.0);
virtual_f_ij.push_back(0.0);
while(true) {
if(virtual_force(i, j) <= 0.0) break;
virtual_initial_integrate(i, j);
if(virtual_force(i, j) <= 0.0) break;
virtual_final_integrate(i, j);
}
{
std::map<int, std::vector<double> >::iterator it = delta_v.find(j);
if (it == delta_v.end()) {
std::vector<double> delta_v0(3,0.0);
delta_v[j] = delta_v0;
}
}
delta_v[i][0] += (virtual_v_i[0] - v[i][0]);
delta_v[i][1] += (virtual_v_i[1] - v[i][1]);
delta_v[i][2] += (virtual_v_i[2] - v[i][2]);
delta_v[j][0] += (virtual_v_j[0] - v[j][0]);
delta_v[j][1] += (virtual_v_j[1] - v[j][1]);
delta_v[j][2] += (virtual_v_j[2] - v[j][2]);
processed_pairs.insert(std::pair<int, int>(i, j));
}
}
}
for (std::map<int, std::vector<double> >::iterator it=delta_v.begin(); it!=delta_v.end(); ++it) {
v[it->first][0] += it->second[0];
v[it->first][1] += it->second[1];
v[it->first][2] += it->second[2];
}
// allocate breakage data
if (breakdata) memory->destroy(breakdata);
memory->create(breakdata,n_break_this_local,10,"FixBreakParticle::breakdata");
// fill breakage data and remove particles
const double eMF = 0.485;
i = ibreak = 0;
while (i < nlocal) {
if (mask[i] & groupbit && flag[i] < 0.0) {
// copy data needed for insertion
vectorCopy3D(x[i],&breakdata[ibreak][0]);
vectorCopy3D(v[i],&breakdata[ibreak][3]);
vectorScalarMult3D(&breakdata[ibreak][3], eMF);
breakdata[ibreak][6] = radius[i];
breakdata[ibreak][7] = -flag[i];
breakdata[ibreak][8] = atom->tag[i];
breakdata[ibreak][9] = (1.0 - eMF*eMF) * 0.5 * rmass[i] * (v[i][0]*v[i][0] + v[i][1]*v[i][1] + v[i][2]*v[i][2]);
++ibreak;
// delete particle
avec->copy(nlocal-1,i,1);
--nlocal;
} else {
++i;
}
}
// update local and global # particles
//NP wait with tags, global map etc
//NP since done by fix insert anyway
atom->nlocal = nlocal;
double rlocal = static_cast<double>(atom->nlocal);
MPI_Allreduce(&rlocal,&atom->natoms,1,MPI_DOUBLE,MPI_SUM,world);
fix_fragments->pre_insert(n_break_this_local, breakdata);
// print stats
print_stats_breakage_during();
}
/* ---------------------------------------------------------------------- */
void FixBreakParticle::print_stats_breakage_during()
{
int step = update->ntimestep;
if (me == 0 && n_break_this > 0) {
if (screen) {
fprintf(screen ,"Particle breakage: broke %d particles (mass %f) at step %d\n - a total of %d particles (mass %f) broken so far.\n",
n_break_this,mass_break_this,step,n_break,mass_break);
}
if (logfile) {
fprintf(logfile,"Particle breakage: broke %d particles (mass %f) at step %d\n - a total of %d particles (mass %f) broken so far.\n",
n_break_this,mass_break_this,step,n_break,mass_break);
}
}
}
/* ---------------------------------------------------------------------- */
int FixBreakParticle::calc_ninsert_this()
{
// number of ptis to insert this timestep
// will effectively insert n_break_this * n_fragments spheres
return n_break_this;
}
/* ----------------------------------------------------------------------
generate new particles at positions where old particles were deleted
function is executed locally on each process as opposed to
FixInsertPack::x_v_omega()
overlap check is not needed since space around broken particles is empty
returns # bodies and # spheres that could actually be inserted
------------------------------------------------------------------------- */
void FixBreakParticle::x_v_omega(int ninsert_this, int &ninserted_this, int &ninserted_spheres_this, double &mass_inserted_this)
{
double pos_ins[3],v_ins[3],omega_ins[3],quat_ins[4];
int nins;
ParticleToInsertFragments *pti;
vectorZeroize3D(omega_ins);
vectorZeroize4D(quat_ins);
// local insertion
double mass_inserted_this_local = 0.;
int ninserted_this_local = 0;
int ninserted_spheres_this_local = 0;
// global insertion
ninserted_this = ninserted_spheres_this = 0;
mass_inserted_this = 0.;
// n_break_this ptis with n_fragments spheres each
// n_break_this * n_fragments spheres to be generated
for (int iparticle = 0; iparticle < n_break_this_local; ++iparticle) {
// get position, velocity and radius of broken particle
vectorCopy3D(&breakdata[iparticle][0],pos_ins);
vectorCopy3D(&breakdata[iparticle][3],v_ins);
// get pti and scale it down with radius of broken particle
pti = static_cast<ParticleToInsertFragments*>(fix_distribution->pti_list[iparticle]);
pti->fix_property_atom_id = fix_collision_factor->id;
nins = pti->set_x_v_omega(pos_ins,v_ins,omega_ins,quat_insert);
// tally stats
ninserted_spheres_this_local += nins;
mass_inserted_this_local += pti->mass_ins;
++ninserted_this_local;
}
// tally stats, have to do this since operation is locally on each process
// as opposed to e.g. FixInsertPack::x_v_omega()
MPI_Sum_Scalar(ninserted_spheres_this_local,ninserted_spheres_this,world);
MPI_Sum_Scalar(ninserted_this_local,ninserted_this,world);
MPI_Sum_Scalar(mass_inserted_this_local,mass_inserted_this,world);
tag_max = atom->tag_max();
}
double FixBreakParticle::virtual_force(int i, int j)
{
virtual_f_ij.clear();
virtual_f_ij.push_back(0.0);
virtual_f_ij.push_back(0.0);
virtual_f_ij.push_back(0.0);
double deltan = 0.0;
if(true) { // hertz/break
const int itype = atom->type[i];
const int jtype = atom->type[j];
double ri = atom->radius[i];
double rj = atom->radius[j];
double reff = ri*rj/(ri+rj);//cdata.is_wall ? cdata.radi : (ri*rj/(ri+rj));
double mi = atom->rmass[i];
double mj = atom->rmass[j];
double meff = mi * mj / (mi + mj);
const double delx = virtual_x_i[0] - virtual_x_j[0];
const double dely = virtual_x_i[1] - virtual_x_j[1];
const double delz = virtual_x_i[2] - virtual_x_j[2];
const double rsq = delx * delx + dely * dely + delz * delz;
const double radsum = ri + rj;
if (rsq < radsum * radsum) {
const double r = sqrt(rsq);
const double rinv = 1.0 / r;
const double enx = delx * rinv;
const double eny = dely * rinv;
const double enz = delz * rinv;
deltan = radsum - r;
double sqrtval = sqrt(reff*deltan);
const int max_type = pair_gran->mpg->max_type();
const double *Y, *nu;
const double * const * e;
Y = static_cast<FixPropertyGlobal*>(modify->find_fix_property("youngsModulus","property/global","peratomtype",max_type,0,style))->get_values();
nu = static_cast<FixPropertyGlobal*>(modify->find_fix_property("poissonsRatio","property/global","peratomtype",max_type,0,style))->get_values();
e = static_cast<FixPropertyGlobal*>(modify->find_fix_property("coefficientRestitution","property/global","peratomtypepair",max_type,max_type,style))->get_array();
const double Yeff = 1./((1.-nu[itype-1]*nu[itype-1])/Y[itype-1]+(1.-nu[jtype-1]*nu[jtype-1])/Y[jtype-1]);
const double loge = log(e[itype-1][jtype-1]);
const double betaeff = loge/sqrt(loge * loge + MY_PI * MY_PI);
const double Sn = 2.*Yeff*sqrtval;
double kn = 4./3. * Yeff * sqrtval;
const double sqrtFiveOverSix = 0.91287092917527685576161630466800355658790782499663875;
double gamman = -2. * sqrtFiveOverSix * betaeff * sqrt(Sn*meff);
// convert Kn and Kt from pressure units to force/distance^2
kn /= force->nktv2p;
// relative translational velocity
const double vr1 = virtual_v_i[0] - virtual_v_j[0];
const double vr2 = virtual_v_i[1] - virtual_v_j[1];
const double vr3 = virtual_v_i[2] - virtual_v_j[2];
const double vn = vr1 * enx + vr2 * eny + vr3 * enz;
const double Fn_damping = -gamman * vn;
const double Fn_contact = kn*deltan;
double Fn = Fn_damping + Fn_contact;
// limit force to avoid the artefact of negative repulsion force
bool limitForce = false;
if (limitForce && (Fn < 0.0) ) {
Fn = 0.0;
}
// apply normal force
if (false) { /*(cdata.is_wall) {
const double Fn_ = Fn * cdata.area_ratio;
i_forces.delta_F[0] = Fn_ * cdata.en[0];
i_forces.delta_F[1] = Fn_ * cdata.en[1];
i_forces.delta_F[2] = Fn_ * cdata.en[2];*/
} else {
virtual_f_ij[0] = Fn * enx;
virtual_f_ij[1] = Fn * eny;
virtual_f_ij[2] = Fn * enz;
}
}
} else { // hooke/break
/// TODO
}
return deltan;
}
void FixBreakParticle::virtual_initial_integrate(int i, int j)
{
const double dtv = update->dt;
const double dtf = 0.5 * update->dt * force->ftm2v;
double *rmass = atom->rmass;
double dtfm;
// update v and x of atoms
dtfm = dtf / rmass[i];
virtual_v_i[0] += dtfm * virtual_f_ij[0];
virtual_v_i[1] += dtfm * virtual_f_ij[1];
virtual_v_i[2] += dtfm * virtual_f_ij[2];
virtual_x_i[0] += dtv * virtual_v_i[0];
virtual_x_i[1] += dtv * virtual_v_i[1];
virtual_x_i[2] += dtv * virtual_v_i[2];
dtfm = dtf / rmass[j];
virtual_v_j[0] -= dtfm * virtual_f_ij[0];
virtual_v_j[1] -= dtfm * virtual_f_ij[1];
virtual_v_j[2] -= dtfm * virtual_f_ij[2];
virtual_x_j[0] += dtv * virtual_v_j[0];
virtual_x_j[1] += dtv * virtual_v_j[1];
virtual_x_j[2] += dtv * virtual_v_j[2];
}
void FixBreakParticle::virtual_final_integrate(int i, int j)
{
const double dtf = 0.5 * update->dt * force->ftm2v;
double *rmass = atom->rmass;
double dtfm;
// update v of atoms
dtfm = dtf / rmass[i];
virtual_v_i[0] += dtfm * virtual_f_ij[0];
virtual_v_i[1] += dtfm * virtual_f_ij[1];
virtual_v_i[2] += dtfm * virtual_f_ij[2];
dtfm = dtf / rmass[j];
virtual_v_j[0] -= dtfm * virtual_f_ij[0];
virtual_v_j[1] -= dtfm * virtual_f_ij[1];
virtual_v_j[2] -= dtfm * virtual_f_ij[2];
}

115
src/fix_break_particle.h Normal file
View File

@ -0,0 +1,115 @@
/* ----------------------------------------------------------------------
LIGGGHTS - LAMMPS Improved for General Granular and Granular Heat
Transfer Simulations
LIGGGHTS is part of the CFDEMproject
www.liggghts.com | www.cfdem.com
Christoph Kloss, christoph.kloss@cfdem.com
Copyright 2009-2012 JKU Linz
Copyright 2012- DCS Computing GmbH, Linz
LIGGGHTS is based on LAMMPS
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
This software is distributed under the GNU General Public License.
See the README file in the top-level directory.
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
Contributing authors:
Christoph Kloss (JKU Linz, DCS Computing GmbH, Linz)
Daniel Queteschiner (JKU Linz)
------------------------------------------------------------------------- */
#ifdef FIX_CLASS
FixStyle(break/particle,FixBreakParticle)
#else
#ifndef LMP_FIX_BREAK_PARTICLE_H
#define LMP_FIX_BREAK_PARTICLE_H
#include "fix_insert.h"
#include <map>
#include <vector>
namespace LAMMPS_NS {
class FixBreakParticle : public FixInsert {
public:
FixBreakParticle(class LAMMPS *, int, char **);
~FixBreakParticle();
void post_create();
void pre_delete(bool unfixflag);
virtual void pre_force(int);
virtual int setmask();
virtual void init();
void init_defaults();
int calc_ninsert_this();
virtual void end_of_step();
protected:
// inherited functions
virtual void calc_insertion_properties();
virtual void pre_insert();
int is_nearby(int);
void x_v_omega(int ninsert_this,int &ninserted_this, int &ninserted_spheres_this, double &mass_inserted_this);
double insertion_fraction();
// functions declared in this class
inline void generate_random(double *pos, double rad_broken,double rad_insert);
void print_stats_breakage_during();
// per breakage flag
class FixPropertyAtom *fix_break;
class FixPropertyAtom *fix_breaker;
class FixPropertyAtom *fix_collision_factor;
// template holding data of the fragments
class FixTemplateFragments *fix_fragments;
// threshold value
double threshold;
int breakage_criterion;
double fMat;
double min_break_rad;
// stats for breakage
int n_break,n_break_this,n_break_this_local;
double mass_break,mass_break_this,mass_break_this_local;
// data of particles to be broken
double **breakdata;
class PairGran *pair_gran;
int dnum;
int deltaMaxOffset;
int siblingOffset;
int collisionFactorOffset;
int impactEnergyOffset;
int tag_max;
double virtual_force(int i, int j);
void virtual_initial_integrate(int i, int j);
void virtual_final_integrate(int i, int j);
std::map<int, std::vector<double> > delta_v;
std::vector<double> virtual_x_i;
std::vector<double> virtual_x_j;
std::vector<double> virtual_v_i;
std::vector<double> virtual_v_j;
std::vector<double> virtual_f_ij;
};
}
#endif
#endif

2
src/fix_insert.cpp
View File

@ -722,6 +722,8 @@ void FixInsert::pre_exchange()
int FixInsert::distribute_ninsert_this(int ninsert_this)
{
if(ninsert_this == 0) return 0;
int me, nprocs, ngap, ninsert_this_local, *ninsert_this_local_all;
double fraction_local, fraction_local_all_sum, *fraction_local_all, *remainder, r, rsum;

299
src/fix_template_fragments.cpp Normal file
View File

@ -0,0 +1,299 @@
/* ----------------------------------------------------------------------
LIGGGHTS - LAMMPS Improved for General Granular and Granular Heat
Transfer Simulations
LIGGGHTS is part of the CFDEMproject
www.liggghts.com | www.cfdem.com
Department for Particule Flow Modelling
Copyright 2014- JKU Linz
LIGGGHTS is based on LAMMPS
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
This software is distributed under the GNU General Public License.
See the README file in the top-level directory.
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
Contributing authors:
Daniel Queteschiner (JKU Linz)
------------------------------------------------------------------------- */
#include "fix_template_fragments.h"
#include "atom.h"
#include "atom_vec.h"
#include "stdlib.h"
#include "vector_liggghts.h"
#include "force.h"
#include "memory.h"
#include "error.h"
#include "particleToInsert_fragments.h"
#include "pair_gran.h"
#include "math_const.h"
#include "particleSizeDistribution.h"
#include "particleSpatialDistribution.h"
using namespace MathConst;
using namespace LAMMPS_NS;
using namespace LMP_PROBABILITY_NS;
/*NL*/#define LMP_DEBUGMODE_FRAGMENTS true
/* ---------------------------------------------------------------------- */
FixTemplateFragments::FixTemplateFragments(LAMMPS *lmp, int narg, char **arg) :
FixTemplateSphere(lmp, narg, arg)
{
if (pdf_density->rand_style() != RANDOM_CONSTANT)
error->all(FLERR,"Fix particletemplate/fragments currently only supports constant density");
if (pdf_radius)
error->fix_error(FLERR,this,"currently does not support keyword 'radius'");
if (strcmp(arg[iarg++],"breakage_index") != 0)
error->fix_error(FLERR,this,"expecting argument 'breakage_index'");
t10_max = atof(arg[iarg++]);
rad_min_pct = atof(arg[iarg++]); // percentage of original size
radiiMassFractions[2] = 0.0;
radiiMassFractions[4] = 0.0;
radiiMassFractions[10] = 0.0;
nspheres = 0;
x_sphere = NULL;
r_sphere = NULL;
}
/* ---------------------------------------------------------------------- */
FixTemplateFragments::~FixTemplateFragments()
{
memory->destroy(x_sphere);
delete [] r_sphere;
}
/* ---------------------------------------------------------------------- */
void FixTemplateFragments::post_create()
{
pair_gran = static_cast<PairGran*>(force->pair_match("gran", 0));
dnum = pair_gran->dnum();
// don't know these in advance ...
volume_expect = 1.0;
mass_expect = volume_expect*expectancy(pdf_density);
r_equiv = 1.0;
r_bound = 0.0;
}
/* ----------------------------------------------------------------------*/
double FixTemplateFragments::min_rad()
{
return 0.0;
}
/* ----------------------------------------------------------------------*/
double FixTemplateFragments::max_rad()
{
return 0.0;
}
/* ----------------------------------------------------------------------*/
double FixTemplateFragments::max_r_bound()
{
return r_bound;
}
/* ----------------------------------------------------------------------*/
// used by pdd to check which template has the most spheres
int FixTemplateFragments::number_spheres()
{
if (LMP_DEBUGMODE_FRAGMENTS) fprintf(screen,"FixTemplateFragments::number_spheres: nspheres = %d\n", nspheres);
return nspheres;
}
/* ----------------------------------------------------------------------*/
void FixTemplateFragments::randomize_single()
{
error->all(FLERR,"re-implement FixTemplateFragments::randomize_single()!");
}
/* ----------------------------------------------------------------------*/
void FixTemplateFragments::init_ptilist(int n_random_max)
{
// n_random_max is local
if (pti_list) error->all(FLERR,"invalid FixTemplateSphere::init_list()");
n_pti_max = n_random_max;
pti_list = (ParticleToInsert**) memory->smalloc(n_pti_max*sizeof(ParticleToInsertFragments*),"pti_list");
memset(pti_list,0,n_pti_max*sizeof(ParticleToInsertFragments*));
}
/* ----------------------------------------------------------------------*/
void FixTemplateFragments::randomize_ptilist(int n_random,int distribution_groupbit)
{
for (int i = 0; i < n_random; ++i) {
nspheres = r_sphere_list[i].size();
delete pti_list[i];
pti_list[i] = new ParticleToInsertFragments(lmp,nspheres);
ParticleToInsertFragments *pti = static_cast<ParticleToInsertFragments*>(pti_list[i]);
pti->density_ins = expectancy(pdf_density);
pti->volume_ins = 0.0;
pti->volume_ins = 0.0;
pti->collision_factor = collision_factor[i];
pti->r_bound_ins = r_bound;
pti->atom_type = atom_type;
for (int j = 0; j < nspheres; ++j) {
pti->radius_ins[j] = r_sphere_list[i][j];
pti->x_ins[j][0] = x_sphere_list[i][j][0];
pti->x_ins[j][1] = x_sphere_list[i][j][1];
pti->x_ins[j][2] = x_sphere_list[i][j][2];
pti->volume_ins += MY_4PI3 * r_sphere_list[i][j] * r_sphere_list[i][j] * r_sphere_list[i][j];
}
pti->mass_ins = pti->volume_ins * pti->density_ins;
vectorZeroize3D(pti->v_ins);
vectorZeroize3D(pti->omega_ins);
pti->groupbit = groupbit | distribution_groupbit; //NP also contains insert_groupbit
}
}
void FixTemplateFragments::pre_insert(int n_total, double **breakdata)
{
// generate particle size/spatial distributions
r_sphere_list.clear();
collision_factor.clear();
r_bound = 0.0;
nspheres = 0;
if (n_total <= 0) {
x_sphere_list.clear();
return;
}
// breakdata[i][6] ... radius
double density = expectancy(pdf_density);
x_sphere_list.resize(n_total);
for (int i = 0; i < n_total; ++i) {
ParticleSizeDistribution psd(breakdata[i][7], density, breakdata[i][6], breakdata[i][6]*rad_min_pct, t10_max);
psd.range_mass_fractions(radiiMassFractions);
std::vector<double> radii;
psd.radii(radiiMassFractions, radii);
r_sphere_list.push_back(radii);
if (r_bound < breakdata[i][6]) {
r_bound = breakdata[i][6];
}
if (nspheres < static_cast<int>(radii.size())) {
nspheres = static_cast<int>(radii.size());
}
// check for overlapping atoms
/// TODO check for overlapping meshes
std::vector<double> ext_radii;
std::vector<std::vector<double> > ext_center;
{
double **x = atom->x;
double *radius = atom->radius;
int nall = atom->nlocal + atom->nghost;
for (int j = 0; j < nall; ++j) { /// TODO optimize: just check neighborlist
const double delx = x[j][0] - breakdata[i][0];
const double dely = x[j][1] - breakdata[i][1];
const double delz = x[j][2] - breakdata[i][2];
const double rsq = delx * delx + dely * dely + delz * delz;
const double radsum = radius[j] + breakdata[i][6];
if (rsq < radsum * radsum) {
std::vector<double> x_j(3, 0.0);
x_j[0] = delx;
x_j[1] = dely;
x_j[2] = delz;
ext_center.push_back(x_j);
ext_radii.push_back(radius[j]);
}
}
}
ParticleSpatialDistribution pxd(random);
x_sphere_list[i].resize(radii.size());
pxd.randomInsertion(breakdata[i][6], radii, x_sphere_list[i], ext_radii, ext_center);
double energy = elastic_energy(r_sphere_list[i], x_sphere_list[i]);
double CF = breakdata[i][9]/energy;
CF = cbrt(CF*CF);
collision_factor.push_back(CF);
}
}
// geometrically correct elastic energy
double FixTemplateFragments::elastic_energy(const std::vector<double> &radius, const std::vector<std::vector<double> > &x)
{
const unsigned int fragments = radius.size();
double energy = 0.0;
int overlaps = 0;
if(fragments > 0) {
const int64_t normalmodel = pair_gran->hashcode() & 0x0f;
const double density = expectancy(pdf_density);
const double *Y, *nu;
const int max_type = pair_gran->mpg->max_type();
const int itype = atom_type - 1;
Y = static_cast<FixPropertyGlobal*>(modify->find_fix_property("youngsModulus","property/global","peratomtype",max_type,0,style))->get_values();
nu = static_cast<FixPropertyGlobal*>(modify->find_fix_property("poissonsRatio","property/global","peratomtype",max_type,0,style))->get_values();
const double Yeff = Y[itype]/( 2. * (1. - nu[itype]*nu[itype]) ); // fragments have same type
for(unsigned int pi = 0; pi < fragments-1; ++pi) {
for(unsigned int pj = pi+1; pj < fragments; ++pj) {
const double delx = x[pj][0] - x[pi][0];
const double dely = x[pj][1] - x[pi][1];
const double delz = x[pj][2] - x[pi][2];
const double rsq = delx * delx + dely * dely + delz * delz;
const double radsum = radius[pi] + radius[pj];
if (rsq < radsum * radsum) {
++overlaps;
const double rij = sqrt(rsq);
const double deltan = radsum - rij;
double reff = radius[pi] * radius[pj] / radsum;
double kn = 0.0;
if (normalmodel == 6) { // hertz/break
kn = 4./3.*Yeff*sqrt(reff * deltan);
energy += 0.4 * kn * deltan * deltan; // 0.4 = (2./5.)
} else if (false) { // hooke/break
const double sqrtval = sqrt(reff);
const double mi = MY_4PI3 * radius[pi] * radius[pi] * radius[pi] * density;
const double mj = MY_4PI3 * radius[pj] * radius[pj] * radius[pj] * density;
const double meff = mi * mj / (mi + mj);
const double charVel = static_cast<FixPropertyGlobal*>(modify->find_fix_property("characteristicVelocity","property/global","scalar",0,0,style))->compute_scalar();
kn = (16./15.) * sqrtval * Yeff * pow(15. * meff * charVel * charVel /(16. * sqrtval * Yeff), 0.2);
energy += 0.5 * kn * deltan * deltan;
}
}
}
}
}
return energy;
}

106
src/fix_template_fragments.h Normal file
View File

@ -0,0 +1,106 @@
/* ----------------------------------------------------------------------
LIGGGHTS - LAMMPS Improved for General Granular and Granular Heat
Transfer Simulations
LIGGGHTS is part of the CFDEMproject
www.liggghts.com | www.cfdem.com
Department for Particule Flow Modelling
Copyright 2014- JKU Linz
LIGGGHTS is based on LAMMPS
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
This software is distributed under the GNU General Public License.
See the README file in the top-level directory.
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
Contributing authors:
Daniel Queteschiner (JKU Linz)
------------------------------------------------------------------------- */
#ifdef FIX_CLASS
FixStyle(particletemplate/fragments,FixTemplateFragments)
#else
#ifndef LMP_FIX_TEMPLATE_FRAGMENTS_H
#define LMP_FIX_TEMPLATE_FRAGMENTS_H
#include <map>
#include <vector>
#include "fix.h"
#include "fix_template_sphere.h"
namespace LAMMPS_NS {
class FixTemplateFragments : public FixTemplateSphere {
public:
FixTemplateFragments(class LAMMPS *, int, char **);
virtual ~FixTemplateFragments();
virtual void post_create();
virtual double min_rad();
virtual double max_rad();
virtual double max_r_bound();
virtual int number_spheres();
// single insertion
virtual void randomize_single();
// multi insertion
virtual void init_ptilist(int);
virtual void randomize_ptilist(int ,int );
void pre_insert(int n_total, double **breakdata);
protected:
// number of spheres in template
int nspheres;
// coords of each sphere with respect to center of mass
double **x_sphere;
// radius of each sphere
double *r_sphere;
// scale factor if read from a file
double scale_fact;
// bounding box
double x_min[3], x_max[3];
// bounding sphere - radius and coordinates with respect to com
double r_bound;
double x_bound[3];
// radius of sphere with equal volume
double r_equiv;
// number of tries for mc
int ntry;
class PairGran *pair_gran;
int dnum;
double rad_min_pct;
double t10_max;
std::map<int, double> radiiMassFractions;
std::vector<std::vector<double> > r_sphere_list;
std::vector<std::vector<std::vector<double> > > x_sphere_list;
std::vector<double> collision_factor;
double elastic_energy(const std::vector<double> &r_sphere, const std::vector<std::vector<double> > &x_sphere);
};
}
#endif
#endif

1
src/math_const.h
View File

@ -22,6 +22,7 @@ namespace MathConst {
static const double MY_2PI = 6.28318530717958647692; // 2pi
static const double MY_3PI = 9.42477796076937971538; // 3pi
static const double MY_4PI = 12.56637061435917295384; // 4pi
static const double MY_4PI3 = 4.18879020478639098462; // 4pi/3
static const double MY_PI2 = 1.57079632679489661923; // pi/2
static const double MY_PI4 = 0.78539816339744830962; // pi/4
static const double MY_PIS = 1.77245385090551602729; // sqrt(pi)

208
src/normal_model_hertz_break.h Normal file
View File

@ -0,0 +1,208 @@
/* ----------------------------------------------------------------------
LIGGGHTS - LAMMPS Improved for General Granular and Granular Heat
Transfer Simulations
LIGGGHTS is part of the CFDEMproject
www.liggghts.com | www.cfdem.com
Christoph Kloss, christoph.kloss@cfdem.com
Copyright 2009-2012 JKU Linz
Copyright 2012- DCS Computing GmbH, Linz
LIGGGHTS is based on LAMMPS
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
This software is distributed under the GNU General Public License.
See the README file in the top-level directory.
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
Contributing authors:
Christoph Kloss (JKU Linz, DCS Computing GmbH, Linz)
Richard Berger (JKU Linz)
Daniel Queteschiner (JKU Linz)
------------------------------------------------------------------------- */
#ifdef NORMAL_MODEL
NORMAL_MODEL(HERTZ_BREAK,hertz/break,6)
#else
#ifndef NORMAL_MODEL_HERTZ_BREAK_H_
#define NORMAL_MODEL_HERTZ_BREAK_H_
#include "contact_models.h"
#include "global_properties.h"
#include "math.h"
namespace LIGGGHTS {
namespace ContactModels
{
template<>
class NormalModel<HERTZ_BREAK> : protected Pointers
{
public:
static const int MASK = CM_REGISTER_SETTINGS | CM_CONNECT_TO_PROPERTIES | CM_COLLISION;
NormalModel(LAMMPS * lmp, IContactHistorySetup * hsetup) : Pointers(lmp),
Yeff(NULL),
Geff(NULL),
betaeff(NULL),
limitForce(false),
displayedSettings(false)
{
history_offset = hsetup->add_history_value("deltaMax", "1");
hsetup->add_history_value("sibling", "1");
hsetup->add_history_value("siblingDeltaMax", "1");
hsetup->add_history_value("collisionFactor", "1");
hsetup->add_history_value("impactEnergy", "1");
/*NL*/ if(comm->me == 0) fprintf(screen, "HERTZ/BREAK loaded\n");
}
void registerSettings(Settings & settings)
{
settings.registerOnOff("tangential_damping", tangential_damping, true);
settings.registerOnOff("limitForce", limitForce);
}
void connectToProperties(PropertyRegistry & registry) {
registry.registerProperty("Yeff", &MODEL_PARAMS::createYeff,"model hertz/break");
registry.registerProperty("Geff", &MODEL_PARAMS::createGeff,"model hertz/break");
registry.registerProperty("betaeff", &MODEL_PARAMS::createBetaEff,"model hertz/break");
registry.connect("Yeff", Yeff,"model hertz/break");
registry.connect("Geff", Geff,"model hertz/break");
registry.connect("betaeff", betaeff,"model hertz/break");
}
// effective exponent for stress-strain relationship
//NP used for area correction of heat transfer
inline double stressStrainExponent()
{
return 1.5;
}
inline void collision(CollisionData & cdata, ForceData & i_forces, ForceData & j_forces)
{
const int itype = cdata.itype;
const int jtype = cdata.jtype;
double ri = cdata.radi;
double rj = cdata.radj;
double reff=cdata.is_wall ? cdata.radi : (ri*rj/(ri+rj));
double meff=cdata.meff;
double deltan = cdata.deltan;
double * const history = &cdata.contact_history[history_offset];
double * const deltaMax = &cdata.contact_history[history_offset];
double sibling = history[1];
double collisionFactor = history[3];
// limit forces between siblings
//NP sibling contact flags must be set after insertion (pre_exchange)
//NP and before this force calculation (pre_force)
if (sibling) {
double * const siblingDeltaMax = &cdata.contact_history[history_offset+2];
if (deltan > *siblingDeltaMax) {
if(*siblingDeltaMax == 0.0) {
*siblingDeltaMax = deltan; // part of the overlap that needs to be scaled
deltan *= collisionFactor;
} else {
deltan -= *siblingDeltaMax;
deltan += *siblingDeltaMax*collisionFactor;
}
} else {
*siblingDeltaMax = deltan;
deltan *= collisionFactor;
}
cdata.shearupdate = false;
cdata.vtr1 = cdata.vtr2 = cdata.vtr3 = 0.0;
}
// detect new contact / maximum overlap
if (*deltaMax == 0.0) { // new contact
*deltaMax = deltan;
history[4] = 0.5 * cdata.vn * cdata.vn; // mass specific impact energy
} else if (*deltaMax > deltan || *deltaMax < 0.0) {
*deltaMax = -deltan;
} else if (*deltaMax <= deltan) {
*deltaMax = deltan;
}
double sqrtval = sqrt(reff*deltan);
double Sn=2.*Yeff[itype][jtype]*sqrtval;
double St=8.*Geff[itype][jtype]*sqrtval;
double kn=4./3.*Yeff[itype][jtype]*sqrtval;
double kt=St;
const double sqrtFiveOverSix = 0.91287092917527685576161630466800355658790782499663875;
double gamman=-2.*sqrtFiveOverSix*betaeff[itype][jtype]*sqrt(Sn*meff);
double gammat=-2.*sqrtFiveOverSix*betaeff[itype][jtype]*sqrt(St*meff);
/*NL*/ //fprintf(screen,"tangential_damping %s\n",tangential_damping?"yes":"no");
if (!tangential_damping) gammat = 0.0;
if (!displayedSettings) {
displayedSettings = true;
/*
if(limitForce)
if(0 == comm->me) fprintf(screen," NormalModel<HERTZ_STIFFNESS>: will limit normal force.\n");
*/
}
// convert Kn and Kt from pressure units to force/distance^2
kn /= force->nktv2p;
kt /= force->nktv2p;
const double Fn_damping = -gamman*cdata.vn;
const double Fn_contact = kn*deltan;
double Fn = Fn_damping + Fn_contact;
//limit force to avoid the artefact of negative repulsion force
if (limitForce && (Fn<0.0) ) {
Fn = 0.0;
}
cdata.Fn = Fn;
cdata.kn = kn;
cdata.kt = kt;
cdata.gamman = gamman;
cdata.gammat = gammat;
// apply normal force
if(cdata.is_wall) {
const double Fn_ = Fn * cdata.area_ratio;
i_forces.delta_F[0] = Fn_ * cdata.en[0];
i_forces.delta_F[1] = Fn_ * cdata.en[1];
i_forces.delta_F[2] = Fn_ * cdata.en[2];
} else {
i_forces.delta_F[0] = cdata.Fn * cdata.en[0];
i_forces.delta_F[1] = cdata.Fn * cdata.en[1];
i_forces.delta_F[2] = cdata.Fn * cdata.en[2];
j_forces.delta_F[0] = -i_forces.delta_F[0];
j_forces.delta_F[1] = -i_forces.delta_F[1];
j_forces.delta_F[2] = -i_forces.delta_F[2];
}
}
void noCollision(ContactData&, ForceData&, ForceData&){}
void beginPass(CollisionData&, ForceData&, ForceData&){}
void endPass(CollisionData&, ForceData&, ForceData&){}
protected:
double ** Yeff;
double ** Geff;
double ** betaeff;
bool tangential_damping;
bool limitForce;
bool displayedSettings;
int history_offset;
};
}
}
#endif
#endif

9
src/pair_gran.h
View File

@ -146,6 +146,15 @@ public:
return offset;
}
int get_history_value_offset(const std::string & name) {
for(std::vector<HistoryArg>::size_type it = 0; it != history_arg.size(); ++it) {
if (history_arg[it].name == name) {
return it;
}
}
return -1;
}
void do_store_contact_forces()
{ store_contact_forces_ = true; }

130
src/particleSizeDistribution.cpp Normal file
View File

@ -0,0 +1,130 @@
/* ----------------------------------------------------------------------
LIGGGHTS - LAMMPS Improved for General Granular and Granular Heat
Transfer Simulations
LIGGGHTS is part of the CFDEMproject
www.liggghts.com | www.cfdem.com
Department for Particule Flow Modelling
Copyright 2014- JKU Linz
LIGGGHTS is based on LAMMPS
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
This software is distributed under the GNU General Public License.
See the README file in the top-level directory.
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
Contributing authors:
Daniel Queteschiner (JKU Linz)
------------------------------------------------------------------------- */
#include <cmath>
#include "math_const.h"
#include "particleSizeDistribution.h"
using namespace LAMMPS_NS;
using namespace MathConst;
ParticleSizeDistribution::ParticleSizeDistribution(double P, double density, double rad_parent, double rad_min, double t10_max) :
breakage_probability(P),
density_(density),
rad_parent_(rad_parent),
rad_min_(rad_min),
t10_max_(t10_max),
t10_(0.0)
{
}
// Note: std::map elements are sorted by their keys.
void ParticleSizeDistribution::mass_fractions(std::map<int, double>& radiiMassFractions)
{
t10_ = t10();
std::map<int,double>::iterator it = radiiMassFractions.begin();
for (; it != radiiMassFractions.end(); ++it) {
it->second = tn(it->first);
}
}
void ParticleSizeDistribution::range_mass_fractions(std::map<int, double>& radiiRangeMassFractions)
{
mass_fractions(radiiRangeMassFractions);
std::map<int,double>::iterator it1 = radiiRangeMassFractions.begin();
std::map<int,double>::iterator it2 = radiiRangeMassFractions.begin();
for (++it2; it2 != radiiRangeMassFractions.end(); ++it1, ++it2) {
it1->second = it1->second - it2->second;
}
}
void ParticleSizeDistribution::radii(const std::map<int, double>& radiiRangeMassFractions, std::vector<double>& radii)
{
const double volume_parent = MY_4PI3 * rad_parent_*rad_parent_*rad_parent_;
const double mass_parent = volume_parent * density_;
double totalMassPool = mass_parent;
std::map<int,double>::const_iterator it1 = radiiRangeMassFractions.begin();
std::map<int,double>::const_iterator it2 = radiiRangeMassFractions.begin();
for (++it2; it1 != radiiRangeMassFractions.end(); ++it1, ++it2) {
double currentParticleRadius = rad_parent_/it1->first;
double currentMassPool = mass_parent * it1->second;
double nextParticleRadiusMax = 0.0;
double nextParticleMassMax = 0.0;
if (it2 != radiiRangeMassFractions.end()) {
nextParticleRadiusMax = rad_parent_/it2->first;
nextParticleMassMax = mass_parent/(it2->first*it2->first*it2->first);
}
// as long as we have enough mass left to form a particle with radius larger than the next section ...
while (currentMassPool > nextParticleMassMax &&
currentParticleRadius > nextParticleRadiusMax &&
currentParticleRadius >= rad_min_) {
double currentParticleMass = MY_4PI3 * currentParticleRadius*currentParticleRadius*currentParticleRadius * density_;
while (currentMassPool >= currentParticleMass) {
radii.push_back(currentParticleRadius);
currentMassPool -= currentParticleMass;
totalMassPool -= currentParticleMass;
}
currentParticleRadius *= 0.999;
}
}
if (totalMassPool > 0.0) {
double radius = cbrt(totalMassPool/(MY_4PI3*density_));
if (radius < rad_min_) {
radius += radii.back();
radii.pop_back();
}
std::vector<double>::iterator it = radii.begin();
while (it != radii.end() && *it > radius)
++it;
radii.insert(it, radius); // slow op; use list?
}
}
double ParticleSizeDistribution::t10()
{
return t10_max_ * breakage_probability;
}
double ParticleSizeDistribution::tn(int n)
{
double alpha = 1.0;
return 1.0 - pow(1.0-t10_, alpha*9.0/(n-1.0));
}

58
src/particleSizeDistribution.h Normal file
View File

@ -0,0 +1,58 @@
/* ----------------------------------------------------------------------
LIGGGHTS - LAMMPS Improved for General Granular and Granular Heat
Transfer Simulations
LIGGGHTS is part of the CFDEMproject
www.liggghts.com | www.cfdem.com
Department for Particule Flow Modelling
Copyright 2014- JKU Linz
LIGGGHTS is based on LAMMPS
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
This software is distributed under the GNU General Public License.
See the README file in the top-level directory.
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
Contributing authors:
Daniel Queteschiner (JKU Linz)
------------------------------------------------------------------------- */
#ifndef LMP_PARTICLE_SIZE_DISTRIBUTION_H
#define LMP_PARTICLE_SIZE_DISTRIBUTION_H
#include <map>
#include <vector>
namespace LAMMPS_NS {
class ParticleSizeDistribution
{
public:
ParticleSizeDistribution(double P, double density, double rad_parent, double rad_min, double t10_max);
void range_mass_fractions(std::map<int, double>& radiiRangeMassFractions);
void radii(const std::map<int, double>& radiiRangeMassFractions, std::vector<double> &radii);
private:
void mass_fractions(std::map<int, double>& radiiMassFractions);
double t10();
double tn(int n);
private:
double breakage_probability;
double density_;
double rad_parent_;
double rad_min_;
double t10_max_;
double t10_;
};
}
#endif // LMP_PARTICLE_SIZE_DISTRIBUTION_H

457
src/particleSpatialDistribution.cpp Normal file
View File

@ -0,0 +1,457 @@
/* ----------------------------------------------------------------------
LIGGGHTS - LAMMPS Improved for General Granular and Granular Heat
Transfer Simulations
LIGGGHTS is part of the CFDEMproject
www.liggghts.com | www.cfdem.com
Department for Particule Flow Modelling
Copyright 2014- JKU Linz
LIGGGHTS is based on LAMMPS
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
This software is distributed under the GNU General Public License.
See the README file in the top-level directory.
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
Contributing authors:
Daniel Queteschiner (JKU Linz)
------------------------------------------------------------------------- */
#include <iostream>
#include <fstream>
#include <sstream>
#include <cmath>
#include <vector>
#include <map>
#include <algorithm>
#include <string>
#include "particleSpatialDistribution.h"
#include "random_park.h"
using namespace LAMMPS_NS;
typedef struct
{
double x;
double y;
double z;
} point;
ParticleSpatialDistribution::ParticleSpatialDistribution(RanPark *rp)
{
RNG = rp;
}
ParticleSpatialDistribution::~ParticleSpatialDistribution()
{
}
bool ParticleSpatialDistribution::isPointInSphere(const std::vector<double> &center, double radius, const std::vector<double> &x, double *dir, double *dist)
{
// dir points from x to center
double delta[3];
delta[0] = center[0] - x[0];
delta[1] = center[1] - x[1];
delta[2] = center[2] - x[2];
double dist_sq = delta[0]*delta[0] + delta[1]*delta[1] + delta[2]*delta[2];
if (dir || dist) {
if (dist) {
double _dist = sqrt(dist_sq);
*dist = _dist;
}
if (dir) {
dir[0] = delta[0];
dir[1] = delta[1];
dir[2] = delta[2];
}
}
if (dist_sq > radius*radius)
return false;
return true;
}
void ParticleSpatialDistribution::randomPointInSphere(double radius, std::vector<double> &x)
{
x[0] = RNG->gaussian();
x[1] = RNG->gaussian();
x[2] = RNG->gaussian();
double length_sq = x[0]*x[0] + x[1]*x[1] + x[2]*x[2];
if (length_sq > 0.0) {
double inverse_length = 1.0/sqrt(length_sq);
double factor = RNG->uniform() * radius * inverse_length;
x[0] *= factor;
x[1] *= factor;
x[2] *= factor;
}
}
#define maxGen 5
#define minRad 0.01
#define clans false
bool create(int generation, int clan, double A, double B, double C, double D, double E, std::multimap<double, point>& pos)
{
// returns false if the radius is less than minRad
// Convert special coords to inversive coords
double rt2o2 = sqrt(2.)/2.;
double a = ( -B + C + D - E) * rt2o2;
double b = ( B - C + D - E) * rt2o2;
double c = ( B + C - D - E) * rt2o2;
double d = (A - B - C - D - E);
double e = ( B + C + D + E) * sqrt(6.)/2.;
// Convert inversive coords to Cartesian coords
// Positive radius required for drawing
double r = 0.05;
if (e > d) {
r = 1./(e - d);
} else {
r = 1./(d - e);
}
/*if(r < minRad) {
return false;
} else*/
{
bool inversionCircle = generation < 0;
bool outerCircle = (generation == 0 && clan == 0);
if (inversionCircle || outerCircle) {
return false;
} else {
point x = {a * r, b * r, c * r};
pos.insert(std::pair<double,point>(r,x));
return true;
}
}
}
double generateItem(int item, double itemVal, int i, double A, double B, double C, double D, double E)
{
switch(i) {
case 0:
if(item == i) return -A;
else return A + itemVal;
case 1:
if(item == i) return -B;
else return B + itemVal;
case 2:
if(item == i) return -C;
else return C + itemVal;
case 3:
if(item == i) return -D;
else return D + itemVal;
default:
if(item == i) return -E;
else return E + itemVal;
}
}
void generate(int generation, int clan, int i, double A, double B, double C, double D, double E, std::multimap<double, point>& pos)
{
double A1 = generateItem(0, A, i, A, B, C, D, E);
double B1 = generateItem(1, B, i, A, B, C, D, E);
double C1 = generateItem(2, C, i, A, B, C, D, E);
double D1 = generateItem(3, D, i, A, B, C, D, E);
double E1 = generateItem(4, E, i, A, B, C, D, E);
// avoid duplicates
if ((i == 0 && (B1 < A1 || C1 < A1 || D1 < A1 || E1 < A1)) ||
(i == 1 && (A1 <=B1 || C1 < B1 || D1 < B1 || E1 < B1)) ||
(i == 2 && (A1 <=C1 || B1 <=C1 || D1 < C1 || E1 < C1)) ||
(i == 3 && (A1 <=D1 || B1 <=D1 || C1 <=D1 || E1 < D1)) ||
(i == 4 && (A1 <=E1 || B1 <=E1 || C1 <=E1 || D1 <=E1))) {
} else {
if (create(generation, clan, A1, B1, C1, D1, E1, pos)) {
if (generation < maxGen) {
for (int j = 0; j < 5; ++j) {
if (j != i) {
generate((generation+1), clan, j, A1, B1, C1, D1, E1, pos);
}
}
}
}
}
}
void ParticleSpatialDistribution::apollonianInsertion(
double radius,
const std::vector<double>& radii,
std::vector<double>& x,
std::vector<double>& y,
std::vector<double>& z)
{
const unsigned int nParticles = radii.size();
if(nParticles == 0) return;
x.clear();
y.clear();
z.clear();
const double r1 = 1.0;
std::multimap<double, point> pos;
// create initial spheres and generate all spheres of their clans
create(0, 0, r1, 0., 0., 0., 0., pos);
create(0, 1, 0., r1, 0., 0., 0., pos);
create(0, 2, 0., 0., r1, 0., 0., pos);
create(0, 3, 0., 0., 0., r1, 0., pos);
create(0, 4, 0., 0., 0., 0., r1, pos);
generate(1, 0, 0, r1, 0., 0., 0., 0., pos);
generate(1, 1, 1, 0., r1, 0., 0., 0., pos);
generate(1, 2, 2, 0., 0., r1, 0., 0., pos);
generate(1, 3, 3, 0., 0., 0., r1, 0., pos);
generate(1, 4, 4, 0., 0., 0., 0., r1, pos);
unsigned int idx = 0;
for (std::multimap<double,point>::reverse_iterator rit = pos.rbegin(); rit != pos.rend(); ++rit, ++idx) {
x.push_back(rit->second.x*radius);
y.push_back(rit->second.y*radius);
z.push_back(rit->second.z*radius);
}
while (x.size() < nParticles) {
x.push_back(0.0);
y.push_back(0.0);
z.push_back(0.0);
}
}
void ParticleSpatialDistribution::randomInsertion(double radius,
const std::vector<double>& radii,
std::vector<std::vector<double> >& x,
const std::vector<double>& ext_radii,
const std::vector<std::vector<double> >& ext_center)
{
const unsigned int nParticles = radii.size();
const unsigned int ext = ext_radii.size();
if (nParticles == 0) return;
for (unsigned int idx = 0; idx < nParticles; ++idx) {
std::vector<double> x_rand(3, 0.0);
bool centerInOccupiedSpace = false;
unsigned int ntry = 0;
do {
randomPointInSphere(radius-radii[idx], x_rand);
++ntry;
// check if point is inside any sphere already inserted -> max overlap < radius
for (unsigned int idx2 = 0; idx2 < idx; ++idx2) {
centerInOccupiedSpace = isPointInSphere(x[idx2], radii[idx2], x_rand);
if (centerInOccupiedSpace)
break;
}
// check if point is inside any external sphere -> max overlap < radius
if (!centerInOccupiedSpace) {
for (unsigned int ext_i = 0; ext_i < ext; ++ext_i) {
std::vector<double> ext_c(3, 0.0);
ext_c[0] = ext_center[ext_i][0];
ext_c[1] = ext_center[ext_i][1];
ext_c[2] = ext_center[ext_i][2];
centerInOccupiedSpace = isPointInSphere(ext_c, ext_radii[ext_i]+radii[idx], x_rand);
if (centerInOccupiedSpace)
break;
}
}
} while (centerInOccupiedSpace && ntry < 200);
x[idx].push_back(x_rand[0]);
x[idx].push_back(x_rand[1]);
x[idx].push_back(x_rand[2]);
}
}
void ParticleSpatialDistribution::relax(
double radius,
const std::vector<double> &radii,
std::vector<std::vector<double> > &x,
const std::vector<double> &ext_radii,
const std::vector<std::vector<double> > &ext_center)
{
const unsigned int nParticles = radii.size();
if (nParticles == 0)
return;
std::vector<double> fx(nParticles, 0.0);
std::vector<double> fy(nParticles, 0.0);
std::vector<double> fz(nParticles, 0.0);
for (unsigned int i = 0; i < 500; ++i) {
calculateForce(radii, x, fx, fy, fz);
applyForce(radius, radii, x, ext_radii, ext_center, fx, fy, fz);
}
}
void ParticleSpatialDistribution::calculateForce(
const std::vector<double> &radii,
const std::vector<std::vector<double> > &x,
std::vector<double> &fx,
std::vector<double> &fy,
std::vector<double> &fz)
{
const unsigned int nParticles = radii.size();
if (nParticles == 0)
return;
fx.reserve(nParticles);
fy.reserve(nParticles);
fz.reserve(nParticles);
std::fill(fx.begin(), fx.end(), 0.0);
std::fill(fy.begin(), fy.end(), 0.0);
std::fill(fz.begin(), fz.end(), 0.0);
for(unsigned int i = 0; i < nParticles-1; ++i) {
for(unsigned int j = i+1; j < nParticles; ++j) {
double delta[3];
delta[0] = x[i][0] - x[j][0];
delta[1] = x[i][1] - x[j][1];
delta[2] = x[i][2] - x[j][2];
double dist_sq = delta[0]*delta[0] + delta[1]*delta[1] + delta[2]*delta[2];
double radsum = radii[i] + radii[j];
if (dist_sq < radsum*radsum) {
const double dist = sqrt(dist_sq);
const double inverse_dist = 1.0 / dist;
const double overlap = radsum - dist;
const double scale = inverse_dist * overlap;
delta[0] *= scale;
delta[1] *= scale;
delta[2] *= scale;
fx[i] += delta[0];
fy[i] += delta[1];
fz[i] += delta[2];
fx[j] -= delta[0];
fy[j] -= delta[1];
fz[j] -= delta[2];
}
}
}
}
void ParticleSpatialDistribution::applyForce(
double radius,
const std::vector<double> &radii,
std::vector<std::vector<double> > &x,
const std::vector<double> &ext_radii,
const std::vector<std::vector<double> > &ext_center,
const std::vector<double> &fx,
const std::vector<double> &fy,
const std::vector<double> &fz)
{
const unsigned int nParticles = radii.size();
const unsigned int ext = ext_radii.size();
if (nParticles == 0)
return;
const double alpha = 0.1;
const std::vector<double> center(3, 0.0);
for (unsigned int i = 0; i < nParticles; ++i) {
x[i][0] += fx[i] * alpha;
x[i][1] += fy[i] * alpha;
x[i][2] += fz[i] * alpha;
double dir[3];
for (unsigned int j = 0; j < ext; ++j) {
std::vector<double> ext_c(3, 0.0);
ext_c[0] = ext_center[j][0];
ext_c[1] = ext_center[j][1];
ext_c[2] = ext_center[j][2];
if (isPointInSphere(ext_c, ext_radii[j] + radii[i], x[i], dir)) {
// moved into contacting sphere
// dir points from pos to center
double dist = sqrt(dir[0]*dir[0] + dir[1]*dir[1] + dir[2]*dir[2]);
double inverse_dist = 1.0/dist;
double scale = -inverse_dist * (ext_radii[j] - dist + radii[i]);
dir[0] *= scale;
dir[1] *= scale;
dir[2] *= scale;
x[i][0] += dir[0];
x[i][1] += dir[1];
x[i][2] += dir[2];
}
}
if (!isPointInSphere(center, radius - radii[i], x[i], dir)) {
// moved outside bounding sphere
double dist = sqrt(dir[0]*dir[0] + dir[1]*dir[1] + dir[2]*dir[2]);
double inverse_dist = 1.0/dist;
double scale = inverse_dist * (dist - (radius - radii[i]));
dir[0] *= scale;
dir[1] *= scale;
dir[2] *= scale;
x[i][0] += dir[0];
x[i][1] += dir[1];
x[i][2] += dir[2];
}
}
}
double ParticleSpatialDistribution::overlap(
const std::vector<double>& radii,
const std::vector<std::vector<double> > &x
) const
{
const unsigned int nParticles = radii.size();
if (nParticles == 0)
return 0.0;
double totalOverlap = 0.0;
for (unsigned int i = 0; i < nParticles-1; ++i) {
for (unsigned int j = i+1; j < nParticles; ++j) {
double delta[3];
delta[0] = x[i][0] - x[j][0];
delta[1] = x[i][1] - x[j][1];
delta[2] = x[i][2] - x[j][2];
double dist_sq = delta[0]*delta[0] + delta[1]*delta[1] + delta[2]*delta[2];
double radsum = radii[i] + radii[j];
if (dist_sq < radsum*radsum)
totalOverlap += radsum - sqrt(dist_sq);
}
}
return totalOverlap;
}

69
src/particleSpatialDistribution.h Normal file
View File

@ -0,0 +1,69 @@
/* ----------------------------------------------------------------------
LIGGGHTS - LAMMPS Improved for General Granular and Granular Heat
Transfer Simulations
LIGGGHTS is part of the CFDEMproject
www.liggghts.com | www.cfdem.com
Department for Particule Flow Modelling
Copyright 2014- JKU Linz
LIGGGHTS is based on LAMMPS
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
This software is distributed under the GNU General Public License.
See the README file in the top-level directory.
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
Contributing authors:
Daniel Queteschiner (JKU Linz)
------------------------------------------------------------------------- */
#ifndef LMP_PARTICLE_SPATIAL_DISTRIBUTION_H
#define LMP_PARTICLE_SPATIAL_DISTRIBUTION_H
namespace LAMMPS_NS {
class RanPark;
class ParticleSpatialDistribution
{
public:
ParticleSpatialDistribution(RanPark *rp);
~ParticleSpatialDistribution();
bool isPointInSphere(const std::vector<double> &center, double radius, const std::vector<double> &x, double *dir=NULL, double *dist=NULL);
void randomPointInSphere(double radius, std::vector<double> &x);
void apollonianInsertion(double radius, const std::vector<double> &radii, std::vector<double> &x, std::vector<double> &y, std::vector<double> &z);
void randomInsertion(double radius,
const std::vector<double> &radii, std::vector<std::vector<double> > &x,
const std::vector<double> &ext_radii, const std::vector<std::vector<double> > &ext_center);
void relax(double radius,
const std::vector<double> &radii, std::vector<std::vector<double> > &x,
const std::vector<double> &ext_radii, const std::vector<std::vector<double> > &ext_center);
void calculateForce(const std::vector<double> &radii,
const std::vector<std::vector<double> > &x,
std::vector<double> &fx, std::vector<double> &fy, std::vector<double> &fz);
void applyForce(double radius,
const std::vector<double> &radii, std::vector<std::vector<double> > &x,
const std::vector<double> &ext_radii, const std::vector<std::vector<double> > &ext_center,
const std::vector<double> &fx, const std::vector<double> &fy, const std::vector<double> &fz);
double overlap(const std::vector<double> &radii, const std::vector<std::vector<double> > &x) const;
private:
RanPark * RNG;
};
}
#endif // SPATIALPARTICLEDISTRIBUTION_H

91
src/particleToInsert_fragments.cpp Normal file
View File

@ -0,0 +1,91 @@
/* ----------------------------------------------------------------------
LIGGGHTS - LAMMPS Improved for General Granular and Granular Heat
Transfer Simulations
LIGGGHTS is part of the CFDEMproject
www.liggghts.com | www.cfdem.com
Department for Particule Flow Modelling
Copyright 2014- JKU Linz
LIGGGHTS is based on LAMMPS
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
This software is distributed under the GNU General Public License.
See the README file in the top-level directory.
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
Contributing authors:
Daniel Queteschiner (JKU Linz)
------------------------------------------------------------------------- */
#include "particleToInsert_fragments.h"
#include "math.h"
#include "math_const.h"
#include "error.h"
#include "update.h"
#include "domain.h"
#include "atom.h"
#include "atom_vec.h"
#include "fix_property_atom.h"
#include "vector_liggghts.h"
#include "modify.h"
using namespace LAMMPS_NS;
using namespace MathConst;
/* ---------------------------------------------------------------------- */
ParticleToInsertFragments::ParticleToInsertFragments(LAMMPS* lmp,int ns) : ParticleToInsert(lmp,ns)
{
collision_factor = 1.0;
}
/* ---------------------------------------------------------------------- */
ParticleToInsertFragments::~ParticleToInsertFragments()
{
}
/* ---------------------------------------------------------------------- */
int ParticleToInsertFragments::insert()
{
// perform the actual insertion
// add particles, set coordinate and radius
// set group mask to "all" plus fix groups
int inserted = 0;
int nfix = modify->nfix;
Fix **fix = modify->fix;
FixPropertyAtom *fix_pa = NULL;
if (!fix_property_atom_id.empty()) {
fix_pa = static_cast<FixPropertyAtom*>(modify->find_fix_id(fix_property_atom_id.c_str()));
}
for(int i = 0; i < nspheres; ++i) {
++inserted;
atom->avec->create_atom(atom_type,x_ins[i]);
int m = atom->nlocal - 1;
atom->mask[m] = 1 | groupbit;
vectorCopy3D(v_ins,atom->v[m]);
vectorCopy3D(omega_ins,atom->omega[m]);
atom->radius[m] = radius_ins[i];
atom->density[m] = density_ins;
atom->rmass[m] = (MY_4PI3) * radius_ins[i]*radius_ins[i]*radius_ins[i] * density_ins;
for (int j = 0; j < nfix; ++j) {
if (fix[j]->create_attribute) fix[j]->set_arrays(m);
}
if (fix_pa) {
fix_pa->vector_atom[m] = collision_factor;
}
}
return inserted;
}

48
src/particleToInsert_fragments.h Normal file
View File

@ -0,0 +1,48 @@
/* ----------------------------------------------------------------------
LIGGGHTS - LAMMPS Improved for General Granular and Granular Heat
Transfer Simulations
LIGGGHTS is part of the CFDEMproject
www.liggghts.com | www.cfdem.com
Department for Particule Flow Modelling
Copyright 2014- JKU Linz
LIGGGHTS is based on LAMMPS
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
This software is distributed under the GNU General Public License.
See the README file in the top-level directory.
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
Contributing authors:
Daniel Queteschiner (JKU Linz)
------------------------------------------------------------------------- */
#ifndef LMP_PARTICLE_TO_INSERT_FRAGMENTS_H
#define LMP_PARTICLE_TO_INSERT_FRAGMENTS_H
#include "particleToInsert.h"
#include <string>
using namespace LAMMPS_NS;
namespace LAMMPS_NS {
class ParticleToInsertFragments : public ParticleToInsert
{
public:
ParticleToInsertFragments(LAMMPS* lmp,int ns = 1);
virtual ~ParticleToInsertFragments();
virtual int insert();
double collision_factor;
std::string fix_property_atom_id;
};
}
#endif