diff --git a/src/fix_break_particle.cpp b/src/fix_break_particle.cpp new file mode 100644 index 00000000..fa82201b --- /dev/null +++ b/src/fix_break_particle.cpp @@ -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 +#include +#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(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(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(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(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(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(n_break_this_local)/static_cast(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(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 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 >::iterator it = delta_v.find(i); + if (it == delta_v.end()) { + std::vector 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 ::iterator, std::multimap::iterator> ret; + ret = processed_pairs.equal_range(j); + for (std::multimap::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 >::iterator it = delta_v.find(j); + if (it == delta_v.end()) { + std::vector 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(i, j)); + } + } + } + + for (std::map >::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(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(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(modify->find_fix_property("youngsModulus","property/global","peratomtype",max_type,0,style))->get_values(); + nu = static_cast(modify->find_fix_property("poissonsRatio","property/global","peratomtype",max_type,0,style))->get_values(); + e = static_cast(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]; +} + + diff --git a/src/fix_break_particle.h b/src/fix_break_particle.h new file mode 100644 index 00000000..70f01ab8 --- /dev/null +++ b/src/fix_break_particle.h @@ -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 +#include + +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 > delta_v; + std::vector virtual_x_i; + std::vector virtual_x_j; + std::vector virtual_v_i; + std::vector virtual_v_j; + std::vector virtual_f_ij; +}; + +} + +#endif +#endif diff --git a/src/fix_insert.cpp b/src/fix_insert.cpp index 0bdc8c37..c94799c4 100644 --- a/src/fix_insert.cpp +++ b/src/fix_insert.cpp @@ -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; diff --git a/src/fix_template_fragments.cpp b/src/fix_template_fragments.cpp new file mode 100644 index 00000000..7165ecb1 --- /dev/null +++ b/src/fix_template_fragments.cpp @@ -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(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(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 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(radii.size())) { + nspheres = static_cast(radii.size()); + } + + // check for overlapping atoms + /// TODO check for overlapping meshes + std::vector ext_radii; + std::vector > 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 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 &radius, const std::vector > &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(modify->find_fix_property("youngsModulus","property/global","peratomtype",max_type,0,style))->get_values(); + nu = static_cast(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(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; +} + diff --git a/src/fix_template_fragments.h b/src/fix_template_fragments.h new file mode 100644 index 00000000..87be46c0 --- /dev/null +++ b/src/fix_template_fragments.h @@ -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 +#include +#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 radiiMassFractions; + std::vector > r_sphere_list; + std::vector > > x_sphere_list; + std::vector collision_factor; + + double elastic_energy(const std::vector &r_sphere, const std::vector > &x_sphere); + +}; + +} + +#endif +#endif diff --git a/src/math_const.h b/src/math_const.h index 89437510..29247904 100644 --- a/src/math_const.h +++ b/src/math_const.h @@ -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) diff --git a/src/normal_model_hertz_break.h b/src/normal_model_hertz_break.h new file mode 100644 index 00000000..41ff7084 --- /dev/null +++ b/src/normal_model_hertz_break.h @@ -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 : 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: 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 diff --git a/src/pair_gran.h b/src/pair_gran.h index d4fa80ef..e7e1aad8 100644 --- a/src/pair_gran.h +++ b/src/pair_gran.h @@ -146,6 +146,15 @@ public: return offset; } + int get_history_value_offset(const std::string & name) { + for(std::vector::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; } diff --git a/src/particleSizeDistribution.cpp b/src/particleSizeDistribution.cpp new file mode 100644 index 00000000..be11dd9d --- /dev/null +++ b/src/particleSizeDistribution.cpp @@ -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 +#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& radiiMassFractions) +{ + t10_ = t10(); + + std::map::iterator it = radiiMassFractions.begin(); + + for (; it != radiiMassFractions.end(); ++it) { + it->second = tn(it->first); + } +} + + +void ParticleSizeDistribution::range_mass_fractions(std::map& radiiRangeMassFractions) +{ + mass_fractions(radiiRangeMassFractions); + + std::map::iterator it1 = radiiRangeMassFractions.begin(); + std::map::iterator it2 = radiiRangeMassFractions.begin(); + + for (++it2; it2 != radiiRangeMassFractions.end(); ++it1, ++it2) { + it1->second = it1->second - it2->second; + } +} + + +void ParticleSizeDistribution::radii(const std::map& radiiRangeMassFractions, std::vector& 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::const_iterator it1 = radiiRangeMassFractions.begin(); + std::map::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::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)); +} diff --git a/src/particleSizeDistribution.h b/src/particleSizeDistribution.h new file mode 100644 index 00000000..c936c2f0 --- /dev/null +++ b/src/particleSizeDistribution.h @@ -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 +#include + +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& radiiRangeMassFractions); + void radii(const std::map& radiiRangeMassFractions, std::vector &radii); + + private: + void mass_fractions(std::map& 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 diff --git a/src/particleSpatialDistribution.cpp b/src/particleSpatialDistribution.cpp new file mode 100644 index 00000000..99935128 --- /dev/null +++ b/src/particleSpatialDistribution.cpp @@ -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 +#include +#include +#include +#include +#include +#include +#include +#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 ¢er, double radius, const std::vector &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 &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& 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(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& 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& radii, + std::vector& x, + std::vector& y, + std::vector& z) +{ + const unsigned int nParticles = radii.size(); + if(nParticles == 0) return; + + x.clear(); + y.clear(); + z.clear(); + + const double r1 = 1.0; + std::multimap 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::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& radii, + std::vector >& x, + const std::vector& ext_radii, + const std::vector >& 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 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 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 &radii, + std::vector > &x, + const std::vector &ext_radii, + const std::vector > &ext_center) +{ + const unsigned int nParticles = radii.size(); + + if (nParticles == 0) + return; + + std::vector fx(nParticles, 0.0); + std::vector fy(nParticles, 0.0); + std::vector 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 &radii, + const std::vector > &x, + std::vector &fx, + std::vector &fy, + std::vector &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 &radii, + std::vector > &x, + const std::vector &ext_radii, + const std::vector > &ext_center, + const std::vector &fx, + const std::vector &fy, + const std::vector &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 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 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& radii, + const std::vector > &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; +} + diff --git a/src/particleSpatialDistribution.h b/src/particleSpatialDistribution.h new file mode 100644 index 00000000..538759bc --- /dev/null +++ b/src/particleSpatialDistribution.h @@ -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 ¢er, double radius, const std::vector &x, double *dir=NULL, double *dist=NULL); + void randomPointInSphere(double radius, std::vector &x); + + void apollonianInsertion(double radius, const std::vector &radii, std::vector &x, std::vector &y, std::vector &z); + + void randomInsertion(double radius, + const std::vector &radii, std::vector > &x, + const std::vector &ext_radii, const std::vector > &ext_center); + + void relax(double radius, + const std::vector &radii, std::vector > &x, + const std::vector &ext_radii, const std::vector > &ext_center); + + void calculateForce(const std::vector &radii, + const std::vector > &x, + std::vector &fx, std::vector &fy, std::vector &fz); + + void applyForce(double radius, + const std::vector &radii, std::vector > &x, + const std::vector &ext_radii, const std::vector > &ext_center, + const std::vector &fx, const std::vector &fy, const std::vector &fz); + + double overlap(const std::vector &radii, const std::vector > &x) const; + + private: + RanPark * RNG; +}; + +} + +#endif // SPATIALPARTICLEDISTRIBUTION_H diff --git a/src/particleToInsert_fragments.cpp b/src/particleToInsert_fragments.cpp new file mode 100644 index 00000000..15a34932 --- /dev/null +++ b/src/particleToInsert_fragments.cpp @@ -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(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; +} diff --git a/src/particleToInsert_fragments.h b/src/particleToInsert_fragments.h new file mode 100644 index 00000000..5647536f --- /dev/null +++ b/src/particleToInsert_fragments.h @@ -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 + +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