Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F81308364
fix_smd_tlsph_reference_configuration.cpp
No One
Temporary
Actions
Download File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Subscribers
None
File Metadata
Details
File Info
Storage
Attached
Created
Thu, Sep 5, 07:02
Size
17 KB
Mime Type
text/x-c
Expires
Sat, Sep 7, 07:02 (2 d)
Engine
blob
Format
Raw Data
Handle
20483134
Attached To
rLAMMPS lammps
fix_smd_tlsph_reference_configuration.cpp
View Options
/* ----------------------------------------------------------------------
*
* *** Smooth Mach Dynamics ***
*
* This file is part of the USER-SMD package for LAMMPS.
* Copyright (2014) Georg C. Ganzenmueller, georg.ganzenmueller@emi.fhg.de
* Fraunhofer Ernst-Mach Institute for High-Speed Dynamics, EMI,
* Eckerstrasse 4, D-79104 Freiburg i.Br, Germany.
*
* This file is based on the FixShearHistory class.
*
* ----------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#include "lattice.h"
#include <mpi.h>
#include <string.h>
#include <stdio.h>
#include "fix_smd_tlsph_reference_configuration.h"
#include "atom.h"
#include "comm.h"
#include "neighbor.h"
#include "neigh_list.h"
#include "force.h"
#include "pair.h"
#include "update.h"
#include "modify.h"
#include "memory.h"
#include "error.h"
#include "domain.h"
#include <Eigen/Eigen>
#include "smd_kernels.h"
#include "smd_math.h"
using namespace Eigen;
using namespace LAMMPS_NS;
using namespace FixConst;
using namespace SMD_Kernels;
using namespace std;
using namespace SMD_Math;
#define DELTA 16384
#define INSERT_PREDEFINED_CRACKS false
/* ---------------------------------------------------------------------- */
FixSMD_TLSPH_ReferenceConfiguration::FixSMD_TLSPH_ReferenceConfiguration(LAMMPS *lmp, int narg, char **arg) :
Fix(lmp, narg, arg) {
if (atom->map_style == 0)
error->all(FLERR, "Pair tlsph with partner list requires an atom map, see atom_modify");
maxpartner = 1;
npartner = NULL;
partner = NULL;
wfd_list = NULL;
wf_list = NULL;
energy_per_bond = NULL;
degradation_ij = NULL;
grow_arrays(atom->nmax);
atom->add_callback(0);
// initialize npartner to 0 so neighbor list creation is OK the 1st time
int nlocal = atom->nlocal;
for (int i = 0; i < nlocal; i++) {
npartner[i] = 0;
}
comm_forward = 14;
updateFlag = 1;
}
/* ---------------------------------------------------------------------- */
FixSMD_TLSPH_ReferenceConfiguration::~FixSMD_TLSPH_ReferenceConfiguration() {
// unregister this fix so atom class doesn't invoke it any more
atom->delete_callback(id, 0);
// delete locally stored arrays
memory->destroy(npartner);
memory->destroy(partner);
memory->destroy(wfd_list);
memory->destroy(wf_list);
memory->destroy(degradation_ij);
memory->destroy(energy_per_bond);
}
/* ---------------------------------------------------------------------- */
int FixSMD_TLSPH_ReferenceConfiguration::setmask() {
int mask = 0;
mask |= PRE_EXCHANGE;
return mask;
}
/* ---------------------------------------------------------------------- */
void FixSMD_TLSPH_ReferenceConfiguration::init() {
if (atom->tag_enable == 0)
error->all(FLERR, "Pair style tlsph requires atoms have IDs");
}
/* ---------------------------------------------------------------------- */
void FixSMD_TLSPH_ReferenceConfiguration::pre_exchange() {
//return;
//printf("in FixSMD_TLSPH_ReferenceConfiguration::pre_exchange()\n");
double **defgrad = atom->smd_data_9;
double *radius = atom->radius;
double *rho = atom->rho;
double *vfrac = atom->vfrac;
double **x = atom->x;
double **x0 = atom->x0;
double *rmass = atom->rmass;
int nlocal = atom->nlocal;
int i, itmp;
int *mask = atom->mask;
if (igroup == atom->firstgroup) {
nlocal = atom->nfirst;
}
int *updateFlag_ptr = (int *) force->pair->extract("smd/tlsph/updateFlag_ptr", itmp);
if (updateFlag_ptr == NULL) {
error->one(FLERR,
"fix FixSMD_TLSPH_ReferenceConfiguration failed to access updateFlag pointer. Check if a pair style exist which calculates this quantity.");
}
int *nn = (int *) force->pair->extract("smd/tlsph/numNeighsRefConfig_ptr", itmp);
if (nn == NULL) {
error->all(FLERR, "FixSMDIntegrateTlsph::updateReferenceConfiguration() failed to access numNeighsRefConfig_ptr array");
}
// sum all update flag across processors
MPI_Allreduce(updateFlag_ptr, &updateFlag, 1, MPI_INT, MPI_MAX, world);
if (updateFlag > 0) {
if (comm->me == 0) {
printf("**** updating ref config at step: %ld\n", update->ntimestep);
}
for (i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
// re-set x0 coordinates
x0[i][0] = x[i][0];
x0[i][1] = x[i][1];
x0[i][2] = x[i][2];
// re-set deformation gradient
defgrad[i][0] = 1.0;
defgrad[i][1] = 0.0;
defgrad[i][2] = 0.0;
defgrad[i][3] = 0.0;
defgrad[i][4] = 1.0;
defgrad[i][5] = 0.0;
defgrad[i][6] = 0.0;
defgrad[i][7] = 0.0;
defgrad[i][8] = 1.0;
/*
* Adjust particle volume as the reference configuration is changed.
* We safeguard against excessive deformations by limiting the adjustment range
* to the intervale J \in [0.9..1.1]
*/
vfrac[i] = rmass[i] / rho[i];
//
if (nn[i] < 15) {
radius[i] *= 1.2;
} // else //{
// radius[i] *= pow(J, 1.0 / domain->dimension);
//}
}
}
// update of reference config could have changed x0, vfrac, radius
// communicate these quantities now to ghosts: x0, vfrac, radius
comm->forward_comm_fix(this);
setup(0);
}
}
/* ----------------------------------------------------------------------
copy partner info from neighbor lists to atom arrays
so can be migrated or stored with atoms
------------------------------------------------------------------------- */
void FixSMD_TLSPH_ReferenceConfiguration::setup(int vflag) {
int i, j, ii, jj, n, inum, jnum;
int *ilist, *jlist, *numneigh, **firstneigh;
double r, h, wf, wfd;
Vector3d dx;
if (updateFlag == 0)
return;
int nlocal = atom->nlocal;
nmax = atom->nmax;
grow_arrays(nmax);
// 1st loop over neighbor list
// calculate npartner for each owned atom
// nlocal_neigh = nlocal when neigh list was built, may be smaller than nlocal
double **x0 = atom->x;
double *radius = atom->radius;
int *mask = atom->mask;
tagint *tag = atom->tag;
NeighList *list = pair->list;
inum = list->inum;
ilist = list->ilist;
numneigh = list->numneigh;
firstneigh = list->firstneigh;
// zero npartner for all current atoms
for (i = 0; i < nlocal; i++)
npartner[i] = 0;
for (ii = 0; ii < inum; ii++) {
i = ilist[ii];
jlist = firstneigh[i];
jnum = numneigh[i];
for (jj = 0; jj < jnum; jj++) {
j = jlist[jj];
j &= NEIGHMASK;
if (INSERT_PREDEFINED_CRACKS) {
if (!crack_exclude(i, j))
continue;
}
dx(0) = x0[i][0] - x0[j][0];
dx(1) = x0[i][1] - x0[j][1];
dx(2) = x0[i][2] - x0[j][2];
r = dx.norm();
h = radius[i] + radius[j];
if (r <= h) {
npartner[i]++;
if (j < nlocal) {
npartner[j]++;
}
}
}
}
maxpartner = 0;
for (i = 0; i < nlocal; i++)
maxpartner = MAX(maxpartner, npartner[i]);
int maxall;
MPI_Allreduce(&maxpartner, &maxall, 1, MPI_INT, MPI_MAX, world);
maxpartner = maxall;
grow_arrays(nmax);
for (i = 0; i < nlocal; i++) {
npartner[i] = 0;
for (jj = 0; jj < maxpartner; jj++) {
wfd_list[i][jj] = 0.0;
wf_list[i][jj] = 0.0;
degradation_ij[i][jj] = 0.0;
energy_per_bond[i][jj] = 0.0;
}
}
for (ii = 0; ii < inum; ii++) {
i = ilist[ii];
jlist = firstneigh[i];
jnum = numneigh[i];
for (jj = 0; jj < jnum; jj++) {
j = jlist[jj];
j &= NEIGHMASK;
dx(0) = x0[i][0] - x0[j][0];
dx(1) = x0[i][1] - x0[j][1];
dx(2) = x0[i][2] - x0[j][2];
r = dx.norm();
h = radius[i] + radius[j];
if (INSERT_PREDEFINED_CRACKS) {
if (!crack_exclude(i, j))
continue;
}
if (r < h) {
spiky_kernel_and_derivative(h, r, domain->dimension, wf, wfd);
partner[i][npartner[i]] = tag[j];
wfd_list[i][npartner[i]] = wfd;
wf_list[i][npartner[i]] = wf;
npartner[i]++;
if (j < nlocal) {
partner[j][npartner[j]] = tag[i];
wfd_list[j][npartner[j]] = wfd;
wf_list[j][npartner[j]] = wf;
npartner[j]++;
}
}
}
}
// count number of particles for which this group is active
// bond statistics
if (update->ntimestep > -1) {
n = 0;
int count = 0;
for (i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
n += npartner[i];
count += 1;
}
}
int nall, countall;
MPI_Allreduce(&n, &nall, 1, MPI_INT, MPI_SUM, world);
MPI_Allreduce(&count, &countall, 1, MPI_INT, MPI_SUM, world);
if (countall < 1) countall = 1;
if (comm->me == 0) {
if (screen) {
printf("\n>>========>>========>>========>>========>>========>>========>>========>>========\n");
fprintf(screen, "TLSPH neighbors:\n");
fprintf(screen, " max # of neighbors for a single particle = %d\n", maxpartner);
fprintf(screen, " average # of neighbors/particle in group tlsph = %g\n", (double) nall / countall);
printf(">>========>>========>>========>>========>>========>>========>>========>>========\n\n");
}
if (logfile) {
fprintf(logfile, "\nTLSPH neighbors:\n");
fprintf(logfile, " max # of neighbors for a single particle = %d\n", maxpartner);
fprintf(logfile, " average # of neighbors/particle in group tlsph = %g\n", (double) nall / countall);
}
}
}
updateFlag = 0; // set update flag to zero after the update
}
/* ----------------------------------------------------------------------
memory usage of local atom-based arrays
------------------------------------------------------------------------- */
double FixSMD_TLSPH_ReferenceConfiguration::memory_usage() {
int nmax = atom->nmax;
int bytes = nmax * sizeof(int);
bytes += nmax * maxpartner * sizeof(tagint); // partner array
bytes += nmax * maxpartner * sizeof(float); // wf_list
bytes += nmax * maxpartner * sizeof(float); // wfd_list
bytes += nmax * maxpartner * sizeof(float); // damage_per_interaction array
bytes += nmax * sizeof(int); // npartner array
return bytes;
}
/* ----------------------------------------------------------------------
allocate local atom-based arrays
------------------------------------------------------------------------- */
void FixSMD_TLSPH_ReferenceConfiguration::grow_arrays(int nmax) {
//printf("in FixSMD_TLSPH_ReferenceConfiguration::grow_arrays\n");
memory->grow(npartner, nmax, "tlsph_refconfig_neigh:npartner");
memory->grow(partner, nmax, maxpartner, "tlsph_refconfig_neigh:partner");
memory->grow(wfd_list, nmax, maxpartner, "tlsph_refconfig_neigh:wfd");
memory->grow(wf_list, nmax, maxpartner, "tlsph_refconfig_neigh:wf");
memory->grow(degradation_ij, nmax, maxpartner, "tlsph_refconfig_neigh:degradation_ij");
memory->grow(energy_per_bond, nmax, maxpartner, "tlsph_refconfig_neigh:damage_onset_strain");
}
/* ----------------------------------------------------------------------
copy values within local atom-based arrays
------------------------------------------------------------------------- */
void FixSMD_TLSPH_ReferenceConfiguration::copy_arrays(int i, int j, int delflag) {
npartner[j] = npartner[i];
for (int m = 0; m < npartner[j]; m++) {
partner[j][m] = partner[i][m];
wfd_list[j][m] = wfd_list[i][m];
wf_list[j][m] = wf_list[i][m];
degradation_ij[j][m] = degradation_ij[i][m];
energy_per_bond[j][m] = energy_per_bond[i][m];
}
}
/* ----------------------------------------------------------------------
pack values in local atom-based arrays for exchange with another proc
------------------------------------------------------------------------- */
int FixSMD_TLSPH_ReferenceConfiguration::pack_exchange(int i, double *buf) {
// NOTE: how do I know comm buf is big enough if extreme # of touching neighs
// Comm::BUFEXTRA may need to be increased
//printf("pack_exchange ...\n");
int m = 0;
buf[m++] = npartner[i];
for (int n = 0; n < npartner[i]; n++) {
buf[m++] = partner[i][n];
buf[m++] = wfd_list[i][n];
buf[m++] = wf_list[i][n];
buf[m++] = degradation_ij[i][n];
buf[m++] = energy_per_bond[i][n];
}
return m;
}
/* ----------------------------------------------------------------------
unpack values in local atom-based arrays from exchange with another proc
------------------------------------------------------------------------- */
int FixSMD_TLSPH_ReferenceConfiguration::unpack_exchange(int nlocal, double *buf) {
if (nlocal == nmax) {
//printf("nlocal=%d, nmax=%d\n", nlocal, nmax);
nmax = nmax / DELTA * DELTA;
nmax += DELTA;
grow_arrays(nmax);
error->message(FLERR,
"in Fixtlsph_refconfigNeighGCG::unpack_exchange: local arrays too small for receiving partner information; growing arrays");
}
//printf("nlocal=%d, nmax=%d\n", nlocal, nmax);
int m = 0;
npartner[nlocal] = static_cast<int>(buf[m++]);
for (int n = 0; n < npartner[nlocal]; n++) {
partner[nlocal][n] = static_cast<tagint>(buf[m++]);
wfd_list[nlocal][n] = static_cast<float>(buf[m++]);
wf_list[nlocal][n] = static_cast<float>(buf[m++]);
degradation_ij[nlocal][n] = static_cast<float>(buf[m++]);
energy_per_bond[nlocal][n] = static_cast<float>(buf[m++]);
}
return m;
}
/* ----------------------------------------------------------------------
pack values in local atom-based arrays for restart file
------------------------------------------------------------------------- */
int FixSMD_TLSPH_ReferenceConfiguration::pack_restart(int i, double *buf) {
int m = 0;
buf[m++] = 4 * npartner[i] + 2;
buf[m++] = npartner[i];
for (int n = 0; n < npartner[i]; n++) {
buf[m++] = partner[i][n];
buf[m++] = wfd_list[i][n];
buf[m++] = wf_list[i][n];
buf[m++] = degradation_ij[i][n];
buf[m++] = energy_per_bond[i][n];
}
return m;
}
/* ----------------------------------------------------------------------
unpack values from atom->extra array to restart the fix
------------------------------------------------------------------------- */
void FixSMD_TLSPH_ReferenceConfiguration::unpack_restart(int nlocal, int nth) {
// ipage = NULL if being called from granular pair style init()
// skip to Nth set of extra values
// double **extra = atom->extra;
//
// int m = 0;
// for (int i = 0; i < nth; i++)
// m += static_cast<int>(extra[nlocal][m]);
// m++;
//
// // allocate new chunks from ipage,dpage for incoming values
//
// npartner[nlocal] = static_cast<int>(extra[nlocal][m++]);
// for (int n = 0; n < npartner[nlocal]; n++) {
// partner[nlocal][n] = static_cast<tagint>(extra[nlocal][m++]);
// }
}
/* ----------------------------------------------------------------------
maxsize of any atom's restart data
------------------------------------------------------------------------- */
int FixSMD_TLSPH_ReferenceConfiguration::maxsize_restart() {
// maxtouch_all = max # of touching partners across all procs
int maxtouch_all;
MPI_Allreduce(&maxpartner, &maxtouch_all, 1, MPI_INT, MPI_MAX, world);
return 4 * maxtouch_all + 2;
}
/* ----------------------------------------------------------------------
size of atom nlocal's restart data
------------------------------------------------------------------------- */
int FixSMD_TLSPH_ReferenceConfiguration::size_restart(int nlocal) {
return 4 * npartner[nlocal] + 2;
}
/* ---------------------------------------------------------------------- */
int FixSMD_TLSPH_ReferenceConfiguration::pack_forward_comm(int n, int *list, double *buf, int pbc_flag, int *pbc) {
int i, j, m;
double *radius = atom->radius;
double *vfrac = atom->vfrac;
double **x0 = atom->x0;
double **defgrad0 = atom->smd_data_9;
//printf("FixSMD_TLSPH_ReferenceConfiguration:::pack_forward_comm\n");
m = 0;
for (i = 0; i < n; i++) {
j = list[i];
buf[m++] = x0[j][0];
buf[m++] = x0[j][1];
buf[m++] = x0[j][2];
buf[m++] = vfrac[j];
buf[m++] = radius[j];
buf[m++] = defgrad0[i][0];
buf[m++] = defgrad0[i][1];
buf[m++] = defgrad0[i][2];
buf[m++] = defgrad0[i][3];
buf[m++] = defgrad0[i][4];
buf[m++] = defgrad0[i][5];
buf[m++] = defgrad0[i][6];
buf[m++] = defgrad0[i][7];
buf[m++] = defgrad0[i][8];
}
return m;
}
/* ---------------------------------------------------------------------- */
void FixSMD_TLSPH_ReferenceConfiguration::unpack_forward_comm(int n, int first, double *buf) {
int i, m, last;
double *radius = atom->radius;
double *vfrac = atom->vfrac;
double **x0 = atom->x0;
double **defgrad0 = atom->smd_data_9;
m = 0;
last = first + n;
for (i = first; i < last; i++) {
x0[i][0] = buf[m++];
x0[i][1] = buf[m++];
x0[i][2] = buf[m++];
vfrac[i] = buf[m++];
radius[i] = buf[m++];
defgrad0[i][0] = buf[m++];
defgrad0[i][1] = buf[m++];
defgrad0[i][2] = buf[m++];
defgrad0[i][3] = buf[m++];
defgrad0[i][4] = buf[m++];
defgrad0[i][5] = buf[m++];
defgrad0[i][6] = buf[m++];
defgrad0[i][7] = buf[m++];
defgrad0[i][8] = buf[m++];
}
}
/* ----------------------------------------------------------------------
routine for excluding bonds across a hardcoded slit crack
Note that everything is scaled by lattice constant l0 to avoid
numerical inaccuracies.
------------------------------------------------------------------------- */
bool FixSMD_TLSPH_ReferenceConfiguration::crack_exclude(int i, int j) {
double **x = atom->x;
double l0 = domain->lattice->xlattice;
// line between pair of atoms i,j
double x1 = x[i][0] / l0;
double y1 = x[i][1] / l0;
double x2 = x[j][0] / l0;
double y2 = x[j][1] / l0;
// hardcoded crack line
double x3 = -0.1 / l0;
double y3 = ((int) 1.0 / l0) + 0.5;
//printf("y3 = %f\n", y3);
double x4 = 0.1 / l0 - 1.0 + 0.1;
double y4 = y3;
bool retVal = DoLineSegmentsIntersect(x1, y1, x2, y2, x3, y3, x4, y4);
return !retVal;
//return 1;
}
Event Timeline
Log In to Comment