/* ----------------------------------------------------------------------
 *
 *                    *** 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.
 *
 * ----------------------------------------------------------------------- */

/* ----------------------------------------------------------------------
 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 "group.h"
#include <math.h>
#include <float.h>
#include <stdlib.h>
#include <string.h>
#include "pair_smd_tlsph.h"
#include "fix_smd_tlsph_reference_configuration.h"
#include "atom.h"
#include "domain.h"
#include "force.h"
#include "update.h"
#include "modify.h"
#include "fix.h"
#include "comm.h"
#include "neighbor.h"
#include "neigh_list.h"
#include "neigh_request.h"
#include "memory.h"
#include "error.h"
#include <stdio.h>
#include <iostream>
#include "math_special.h"
#include <map>
#include "update.h"
#include <Eigen/Eigen>
#include "smd_material_models.h"
#include "smd_kernels.h"
#include "smd_math.h"
using namespace SMD_Kernels;
using namespace Eigen;
using namespace std;
using namespace LAMMPS_NS;
using namespace SMD_Math;

#define JAUMANN false
#define DETF_MIN 0.2 // maximum compression deformation allow
#define DETF_MAX 2.0 // maximum tension deformation allowed
#define TLSPH_DEBUG 0
#define PLASTIC_STRAIN_AVERAGE_WINDOW 100.0

/* ---------------------------------------------------------------------- */

PairTlsph::PairTlsph(LAMMPS *lmp) :
		Pair(lmp) {

	onerad_dynamic = onerad_frozen = maxrad_dynamic = maxrad_frozen = NULL;

	failureModel = NULL;
	strengthModel = eos = NULL;

	nmax = 0; // make sure no atom on this proc such that initial memory allocation is correct
	Fdot = Fincr = K = PK1 = NULL;
	R = FincrInv = W = D = NULL;
	detF = NULL;
	smoothVelDifference = NULL;
	numNeighsRefConfig = NULL;
	CauchyStress = NULL;
	hourglass_error = NULL;
	Lookup = NULL;
	particle_dt = NULL;

	updateFlag = 0;
	first = true;
	dtCFL = 0.0; // initialize dtCFL so it is set to safe value if extracted on zero-th timestep

	comm_forward = 22; // this pair style communicates 20 doubles to ghost atoms : PK1 tensor + F tensor + shepardWeight
	fix_tlsph_reference_configuration = NULL;

	cut_comm = MAX(neighbor->cutneighmax, comm->cutghostuser); // cutoff radius within which ghost atoms are communicated.
}

/* ---------------------------------------------------------------------- */

PairTlsph::~PairTlsph() {
	//printf("in PairTlsph::~PairTlsph()\n");

	if (allocated) {
		memory->destroy(setflag);
		memory->destroy(cutsq);
		memory->destroy(strengthModel);
		memory->destroy(eos);
		memory->destroy(Lookup);

		delete[] onerad_dynamic;
		delete[] onerad_frozen;
		delete[] maxrad_dynamic;
		delete[] maxrad_frozen;

		delete[] Fdot;
		delete[] Fincr;
		delete[] K;
		delete[] detF;
		delete[] PK1;
		delete[] smoothVelDifference;
		delete[] R;
		delete[] FincrInv;
		delete[] W;
		delete[] D;
		delete[] numNeighsRefConfig;
		delete[] CauchyStress;
		delete[] hourglass_error;
		delete[] particle_dt;

		delete[] failureModel;
	}
}

/* ----------------------------------------------------------------------
 *
 * use half neighbor list to re-compute shape matrix
 *
 ---------------------------------------------------------------------- */

void PairTlsph::PreCompute() {
	tagint *mol = atom->molecule;
	double *vfrac = atom->vfrac;
	double *radius = atom->radius;
	double **x0 = atom->x0;
	double **x = atom->x;
	double **v = atom->vest; // extrapolated velocities corresponding to current positions
	double **vint = atom->v; // Velocity-Verlet algorithm velocities
	double *damage = atom->damage;
	tagint *tag = atom->tag;
	int *type = atom->type;
	int nlocal = atom->nlocal;
	int jnum, jj, i, j, itype, idim;

	tagint **partner = ((FixSMD_TLSPH_ReferenceConfiguration *) modify->fix[ifix_tlsph])->partner;
	int *npartner = ((FixSMD_TLSPH_ReferenceConfiguration *) modify->fix[ifix_tlsph])->npartner;
	float **wfd_list = ((FixSMD_TLSPH_ReferenceConfiguration *) modify->fix[ifix_tlsph])->wfd_list;
	float **wf_list = ((FixSMD_TLSPH_ReferenceConfiguration *) modify->fix[ifix_tlsph])->wf_list;
	float **degradation_ij = ((FixSMD_TLSPH_ReferenceConfiguration *) modify->fix[ifix_tlsph])->degradation_ij;
	double r0, r0Sq, wf, wfd, h, irad, voli, volj, scale, shepardWeight;
	Vector3d dx, dx0, dv, g;
	Matrix3d Ktmp, Ftmp, Fdottmp, L, U, eye;
	Vector3d vi, vj, vinti, vintj, xi, xj, x0i, x0j, dvint;
	int periodic = (domain->xperiodic || domain->yperiodic || domain->zperiodic);
	bool status;
	Matrix3d F0;

	eye.setIdentity();

	for (i = 0; i < nlocal; i++) {

		itype = type[i];
		if (setflag[itype][itype] == 1) {

			K[i].setZero();
			Fincr[i].setZero();
			Fdot[i].setZero();
			numNeighsRefConfig[i] = 0;
			smoothVelDifference[i].setZero();
			hourglass_error[i] = 0.0;

			if (mol[i] < 0) { // valid SPH particle have mol > 0
				continue;
			}

			// initialize aveage mass density
			h = 2.0 * radius[i];
			r0 = 0.0;
			spiky_kernel_and_derivative(h, r0, domain->dimension, wf, wfd);
			shepardWeight = wf * voli;

			jnum = npartner[i];
			irad = radius[i];
			voli = vfrac[i];

			// initialize Eigen data structures from LAMMPS data structures
			for (idim = 0; idim < 3; idim++) {
				xi(idim) = x[i][idim];
				x0i(idim) = x0[i][idim];
				vi(idim) = v[i][idim];
				vinti(idim) = vint[i][idim];
			}

			for (jj = 0; jj < jnum; jj++) {

				if (partner[i][jj] == 0)
					continue;
				j = atom->map(partner[i][jj]);
				if (j < 0) { //			// check if lost a partner without first breaking bond
					partner[i][jj] = 0;
					continue;
				}

				if (mol[j] < 0) { // particle has failed. do not include it for computing any property
					continue;
				}

				if (mol[i] != mol[j]) {
					continue;
				}

				// initialize Eigen data structures from LAMMPS data structures
				for (idim = 0; idim < 3; idim++) {
					xj(idim) = x[j][idim];
					x0j(idim) = x0[j][idim];
					vj(idim) = v[j][idim];
					vintj(idim) = vint[j][idim];
				}
				dx0 = x0j - x0i;
				dx = xj - xi;

				if (periodic)
					domain->minimum_image(dx0(0), dx0(1), dx0(2));

				r0Sq = dx0.squaredNorm();
				h = irad + radius[j];

				r0 = sqrt(r0Sq);
				volj = vfrac[j];

				// distance vectors in current and reference configuration, velocity difference
				dv = vj - vi;
				dvint = vintj - vinti;

				// scale the interaction according to the damage variable
				scale = 1.0 - degradation_ij[i][jj];
				wf = wf_list[i][jj] * scale;
				wfd = wfd_list[i][jj] * scale;
				g = (wfd / r0) * dx0;

				/* build matrices */
				Ktmp = -g * dx0.transpose();
				Fdottmp = -dv * g.transpose();
				Ftmp = -(dx - dx0) * g.transpose();

				K[i] += volj * Ktmp;
				Fdot[i] += volj * Fdottmp;
				Fincr[i] += volj * Ftmp;
				shepardWeight += volj * wf;
				smoothVelDifference[i] += volj * wf * dvint;
				numNeighsRefConfig[i]++;
			} // end loop over j

			// normalize average velocity field around an integration point
			if (shepardWeight > 0.0) {
				smoothVelDifference[i] /= shepardWeight;
			} else {
				smoothVelDifference[i].setZero();
			}

			pseudo_inverse_SVD(K[i]);
			Fdot[i] *= K[i];
			Fincr[i] *= K[i];
			Fincr[i] += eye;

			if (JAUMANN) {
				R[i].setIdentity(); // for Jaumann stress rate, we do not need a subsequent rotation back into the reference configuration
			} else {
				status = PolDec(Fincr[i], R[i], U, false); // polar decomposition of the deformation gradient, F = R * U
				if (!status) {
					error->message(FLERR, "Polar decomposition of deformation gradient failed.\n");
					mol[i] = -1;
				} else {
					Fincr[i] = R[i] * U;
				}
			}

			detF[i] = Fincr[i].determinant();
			FincrInv[i] = Fincr[i].inverse();

			// velocity gradient
			L = Fdot[i] * FincrInv[i];

			// symmetric (D) and asymmetric (W) parts of L
			D[i] = 0.5 * (L + L.transpose());
			W[i] = 0.5 * (L - L.transpose()); // spin tensor:: need this for Jaumann rate

			// unrotated rate-of-deformation tensor d, see right side of Pronto2d, eqn.(2.1.7)
			// convention: unrotated frame is that one, where the true rotation of an integration point has been subtracted.
			// stress in the unrotated frame of reference is denoted sigma (stress seen by an observer doing rigid body rotations along with the material)
			// stress in the true frame of reference (a stationary observer) is denoted by T, "true stress"
			D[i] = (R[i].transpose() * D[i] * R[i]).eval();

			// limit strain rate
			//double limit = 1.0e-3 * Lookup[SIGNAL_VELOCITY][itype] / radius[i];
			//D[i] = LimitEigenvalues(D[i], limit);

			/*
			 * make sure F stays within some limits
			 */

			if ((detF[i] < DETF_MIN) || (detF[i] > DETF_MAX) || (numNeighsRefConfig[i] == 0)) {
				printf("deleting particle [%d] because det(F)=%f is outside stable range %f -- %f \n", tag[i],
						Fincr[i].determinant(),
						DETF_MIN, DETF_MAX);
				printf("nn = %d, damage=%f\n", numNeighsRefConfig[i], damage[i]);
				cout << "Here is matrix F:" << endl << Fincr[i] << endl;
				cout << "Here is matrix F-1:" << endl << FincrInv[i] << endl;
				cout << "Here is matrix K-1:" << endl << K[i] << endl;
				cout << "Here is matrix K:" << endl << K[i].inverse() << endl;
				cout << "Here is det of K" << endl << (K[i].inverse()).determinant() << endl;
				cout << "Here is matrix R:" << endl << R[i] << endl;
				cout << "Here is det of R" << endl << R[i].determinant() << endl;
				cout << "Here is matrix U:" << endl << U << endl;
				mol[i] = -1;
				//error->one(FLERR, "");
			}

			if (mol[i] < 0) {
				D[i].setZero();
				Fdot[i].setZero();
				Fincr[i].setIdentity();
				smoothVelDifference[i].setZero();
				detF[i] = 1.0;
				K[i].setIdentity();

				vint[i][0] = 0.0;
				vint[i][1] = 0.0;
				vint[i][2] = 0.0;
			}
		} // end loop over i
	} // end check setflag
}

/* ---------------------------------------------------------------------- */

void PairTlsph::compute(int eflag, int vflag) {

	if (atom->nmax > nmax) {
		nmax = atom->nmax;
		delete[] Fdot;
		Fdot = new Matrix3d[nmax]; // memory usage: 9 doubles
		delete[] Fincr;
		Fincr = new Matrix3d[nmax]; // memory usage: 9 doubles
		delete[] K;
		K = new Matrix3d[nmax]; // memory usage: 9 doubles
		delete[] PK1;
		PK1 = new Matrix3d[nmax]; // memory usage: 9 doubles; total 5*9=45 doubles
		delete[] detF;
		detF = new double[nmax]; // memory usage: 1 double; total 46 doubles
		delete[] smoothVelDifference;
		smoothVelDifference = new Vector3d[nmax]; // memory usage: 3 doubles; total 49 doubles
		delete[] R;
		R = new Matrix3d[nmax]; // memory usage: 9 doubles; total 67 doubles
		delete[] FincrInv;
		FincrInv = new Matrix3d[nmax]; // memory usage: 9 doubles; total 85 doubles
		delete[] W;
		W = new Matrix3d[nmax]; // memory usage: 9 doubles; total 94 doubles
		delete[] D;
		D = new Matrix3d[nmax]; // memory usage: 9 doubles; total 103 doubles
		delete[] numNeighsRefConfig;
		numNeighsRefConfig = new int[nmax]; // memory usage: 1 int; total 108 doubles
		delete[] CauchyStress;
		CauchyStress = new Matrix3d[nmax]; // memory usage: 9 doubles; total 118 doubles
		delete[] hourglass_error;
		hourglass_error = new double[nmax];
		delete[] particle_dt;
		particle_dt = new double[nmax];
	}

	if (first) { // return on first call, because reference connectivity lists still needs to be built. Also zero quantities which are otherwise undefined.
		first = false;

		for (int i = 0; i < atom->nlocal; i++) {
			Fincr[i].setZero();
			detF[i] = 0.0;
			smoothVelDifference[i].setZero();
			D[i].setZero();
			numNeighsRefConfig[i] = 0;
			CauchyStress[i].setZero();
			hourglass_error[i] = 0.0;
			particle_dt[i] = 0.0;
		}

		return;
	}

	/*
	 * calculate deformations and rate-of-deformations
	 */
	PairTlsph::PreCompute();

	/*
	 * calculate stresses from constitutive models
	 */
	PairTlsph::AssembleStress();

	/*
	 * QUANTITIES ABOVE HAVE ONLY BEEN CALCULATED FOR NLOCAL PARTICLES.
	 * NEED TO DO A FORWARD COMMUNICATION TO GHOST ATOMS NOW
	 */
	comm->forward_comm_pair(this);

	/*
	 * compute forces between particles
	 */
	updateFlag = 0;
	ComputeForces(eflag, vflag);
}

void PairTlsph::ComputeForces(int eflag, int vflag) {
	tagint *mol = atom->molecule;
	double **x = atom->x;
	double **v = atom->vest;
	double **x0 = atom->x0;
	double **f = atom->f;
	double *vfrac = atom->vfrac;
	double *de = atom->de;
	double *rmass = atom->rmass;
	double *radius = atom->radius;
	double *damage = atom->damage;
	double *plastic_strain = atom->eff_plastic_strain;
	int *type = atom->type;
	int nlocal = atom->nlocal;
	int i, j, jj, jnum, itype, idim;
	double r, hg_mag, wf, wfd, h, r0, r0Sq, voli, volj;
	double delVdotDelR, visc_magnitude, deltaE, mu_ij, hg_err, gamma_dot_dx, delta, scale;
	double strain1d, strain1d_max, softening_strain, shepardWeight;
	char str[128];
	Vector3d fi, fj, dx0, dx, dv, f_stress, f_hg, dxp_i, dxp_j, gamma, g, gamma_i, gamma_j, x0i, x0j;
	Vector3d xi, xj, vi, vj, f_visc, sumForces, f_spring;
	int periodic = (domain->xperiodic || domain->yperiodic || domain->zperiodic);

	tagint **partner = ((FixSMD_TLSPH_ReferenceConfiguration *) modify->fix[ifix_tlsph])->partner;
	int *npartner = ((FixSMD_TLSPH_ReferenceConfiguration *) modify->fix[ifix_tlsph])->npartner;
	float **wfd_list = ((FixSMD_TLSPH_ReferenceConfiguration *) modify->fix[ifix_tlsph])->wfd_list;
	float **wf_list = ((FixSMD_TLSPH_ReferenceConfiguration *) modify->fix[ifix_tlsph])->wf_list;
	float **degradation_ij = ((FixSMD_TLSPH_ReferenceConfiguration *) modify->fix[ifix_tlsph])->degradation_ij;
	float **energy_per_bond = ((FixSMD_TLSPH_ReferenceConfiguration *) modify->fix[ifix_tlsph])->energy_per_bond;
	Matrix3d eye;
	eye.setIdentity();

	if (eflag || vflag)
		ev_setup(eflag, vflag);
	else
		evflag = vflag_fdotr = 0;

	/*
	 * iterate over pairs of particles i, j and assign forces using PK1 stress tensor
	 */

	//updateFlag = 0;
	hMin = 1.0e22;
	dtRelative = 1.0e22;

	for (i = 0; i < nlocal; i++) {

		if (mol[i] < 0) {
			continue; // Particle i is not a valid SPH particle (anymore). Skip all interactions with this particle.
		}

		itype = type[i];
		jnum = npartner[i];
		voli = vfrac[i];

		// initialize aveage mass density
		h = 2.0 * radius[i];
		r = 0.0;
		spiky_kernel_and_derivative(h, r, domain->dimension, wf, wfd);
		shepardWeight = wf * voli;

		for (idim = 0; idim < 3; idim++) {
			x0i(idim) = x0[i][idim];
			xi(idim) = x[i][idim];
			vi(idim) = v[i][idim];
		}

		for (jj = 0; jj < jnum; jj++) {
			if (partner[i][jj] == 0)
				continue;
			j = atom->map(partner[i][jj]);
			if (j < 0) { //			// check if lost a partner without first breaking bond
				partner[i][jj] = 0;
				continue;
			}

			if (mol[j] < 0) {
				continue; // Particle j is not a valid SPH particle (anymore). Skip all interactions with this particle.
			}

			if (mol[i] != mol[j]) {
				continue;
			}

			if (type[j] != itype) {
				sprintf(str, "particle pair is not of same type!");
				error->all(FLERR, str);
			}

			for (idim = 0; idim < 3; idim++) {
				x0j(idim) = x0[j][idim];
				xj(idim) = x[j][idim];
				vj(idim) = v[j][idim];
			}

			if (periodic)
				domain->minimum_image(dx0(0), dx0(1), dx0(2));

			// check that distance between i and j (in the reference config) is less than cutoff
			dx0 = x0j - x0i;
			r0Sq = dx0.squaredNorm();
			h = radius[i] + radius[j];
			hMin = MIN(hMin, h);
			r0 = sqrt(r0Sq);
			volj = vfrac[j];

			// distance vectors in current and reference configuration, velocity difference
			dx = xj - xi;
			dv = vj - vi;
			r = dx.norm(); // current distance

			// scale the interaction according to the damage variable
			scale = 1.0 - degradation_ij[i][jj];
			wf = wf_list[i][jj] * scale;
			wfd = wfd_list[i][jj] * scale;

			g = (wfd / r0) * dx0; // uncorrected kernel gradient

			/*
			 * force contribution -- note that the kernel gradient correction has been absorbed into PK1
			 */

			f_stress = -voli * volj * (PK1[i] + PK1[j]) * g;

			/*
			 * artificial viscosity
			 */
			delVdotDelR = dx.dot(dv) / (r + 0.1 * h); // project relative velocity onto unit particle distance vector [m/s]
			LimitDoubleMagnitude(delVdotDelR, 0.01 * Lookup[SIGNAL_VELOCITY][itype]);
			mu_ij = h * delVdotDelR / (r + 0.1 * h); // units: [m * m/s / m = m/s]
			visc_magnitude = (-Lookup[VISCOSITY_Q1][itype] * Lookup[SIGNAL_VELOCITY][itype] * mu_ij
					+ Lookup[VISCOSITY_Q2][itype] * mu_ij * mu_ij) / Lookup[REFERENCE_DENSITY][itype]; // units: m^5/(s^2 kg))
			f_visc = rmass[i] * rmass[j] * visc_magnitude * wfd * dx / (r + 1.0e-2 * h); // units: kg^2 * m^5/(s^2 kg) * m^-4 = kg m / s^2 = N

			/*
			 * hourglass deviation of particles i and j
			 */

			gamma = 0.5 * (Fincr[i] + Fincr[j]) * dx0 - dx;
			hg_err = gamma.norm() / r0;
			hourglass_error[i] += volj * wf * hg_err;

			/* SPH-like hourglass formulation */

			if (MAX(plastic_strain[i], plastic_strain[j]) > 1.0e-3) {
				/*
				 * viscous hourglass formulation for particles with plastic deformation
				 */
				delta = gamma.dot(dx);
				if (delVdotDelR * delta < 0.0) {
					hg_err = MAX(hg_err, 0.05); // limit hg_err to avoid numerical instabilities
					hg_mag = -hg_err * Lookup[HOURGLASS_CONTROL_AMPLITUDE][itype] * Lookup[SIGNAL_VELOCITY][itype] * mu_ij
							/ Lookup[REFERENCE_DENSITY][itype]; // this has units of pressure
				} else {
					hg_mag = 0.0;
				}
				f_hg = rmass[i] * rmass[j] * hg_mag * wfd * dx / (r + 1.0e-2 * h);

			} else {
				/*
				 * stiffness hourglass formulation for particle in the elastic regime
				 */

				gamma_dot_dx = gamma.dot(dx); // project hourglass error vector onto pair distance vector
				LimitDoubleMagnitude(gamma_dot_dx, 0.1 * r); // limit projected vector to avoid numerical instabilities
				delta = 0.5 * gamma_dot_dx / (r + 0.1 * h); // delta has dimensions of [m]
				hg_mag = Lookup[HOURGLASS_CONTROL_AMPLITUDE][itype] * delta / (r0Sq + 0.01 * h * h); // hg_mag has dimensions [m^(-1)]
				hg_mag *= -voli * volj * wf * Lookup[YOUNGS_MODULUS][itype]; // hg_mag has dimensions [J*m^(-1)] = [N]
				f_hg = (hg_mag / (r + 0.01 * h)) * dx;
			}

			// scale hourglass force with damage
			f_hg *= (1.0 - damage[i]) * (1.0 - damage[j]);

			// sum stress, viscous, and hourglass forces
			sumForces = f_stress + f_visc + f_hg; // + f_spring;

			// energy rate -- project velocity onto force vector
			deltaE = 0.5 * sumForces.dot(dv);

			// apply forces to pair of particles
			f[i][0] += sumForces(0);
			f[i][1] += sumForces(1);
			f[i][2] += sumForces(2);
			de[i] += deltaE;

			// tally atomistic stress tensor
			if (evflag) {
				ev_tally_xyz(i, j, nlocal, 0, 0.0, 0.0, sumForces(0), sumForces(1), sumForces(2), dx(0), dx(1), dx(2));
			}

			shepardWeight += wf * volj;

			// check if a particle has moved too much w.r.t another particle
			if (r > r0) {
				if (update_method == UPDATE_CONSTANT_THRESHOLD) {
					if (r - r0 > update_threshold) {
						updateFlag = 1;
					}
				} else if (update_method == UPDATE_PAIRWISE_RATIO) {
					if ((r - r0) / h > update_threshold) {
						updateFlag = 1;
					}
				}
			}

			if (failureModel[itype].failure_max_pairwise_strain) {

				strain1d = (r - r0) / r0;
				strain1d_max = Lookup[FAILURE_MAX_PAIRWISE_STRAIN_THRESHOLD][itype];
				softening_strain = 2.0 * strain1d_max;

				if (strain1d > strain1d_max) {
					degradation_ij[i][jj] = (strain1d - strain1d_max) / softening_strain;
				} else {
					degradation_ij[i][jj] = 0.0;
				}

				if (degradation_ij[i][jj] >= 1.0) { // delete interaction if fully damaged
					partner[i][jj] = 0;
				}
			}

			if (failureModel[itype].failure_energy_release_rate) {

				// integration approach
				energy_per_bond[i][jj] += update->dt * f_stress.dot(dv) / (voli * volj);
				double Vic = (2.0 / 3.0) * h * h * h; // interaction volume for 2d plane strain
				double critical_energy_per_bond = Lookup[CRITICAL_ENERGY_RELEASE_RATE][itype] / (2.0 * Vic);

				if (energy_per_bond[i][jj] > critical_energy_per_bond) {
					//degradation_ij[i][jj] = 1.0;
					partner[i][jj] = 0;
				}
			}

			if (failureModel[itype].integration_point_wise) {

				strain1d = (r - r0) / r0;

				if (strain1d > 0.0) {

					if ((damage[i] == 1.0) && (damage[j] == 1.0)) {
						// check if damage_onset is already defined
						if (energy_per_bond[i][jj] == 0.0) { // pair damage not defined yet
							energy_per_bond[i][jj] = strain1d;
						} else { // damage initiation strain already defined
							strain1d_max = energy_per_bond[i][jj];
							softening_strain = 2.0 * strain1d_max;

							if (strain1d > strain1d_max) {
								degradation_ij[i][jj] = (strain1d - strain1d_max) / softening_strain;
							} else {
								degradation_ij[i][jj] = 0.0;
							}
						}
					}

					if (degradation_ij[i][jj] >= 1.0) { // delete interaction if fully damaged
						partner[i][jj] = 0;
					}

				} else {
					degradation_ij[i][jj] = 0.0;
				} // end failureModel[itype].integration_point_wise

			}

		} // end loop over jj neighbors of i

		if (shepardWeight != 0.0) {
			hourglass_error[i] /= shepardWeight;
		}

	} // end loop over i

	if (vflag_fdotr)
		virial_fdotr_compute();
}

/* ----------------------------------------------------------------------
 assemble unrotated stress tensor using deviatoric and pressure components.
 Convert to corotational Cauchy stress, then to PK1 stress and apply
 shape matrix correction
 ------------------------------------------------------------------------- */
void PairTlsph::AssembleStress() {
	tagint *mol = atom->molecule;
	double *eff_plastic_strain = atom->eff_plastic_strain;
	double *eff_plastic_strain_rate = atom->eff_plastic_strain_rate;
	double **tlsph_stress = atom->smd_stress;
	int *type = atom->type;
	double *radius = atom->radius;
	double *damage = atom->damage;
	double *rmass = atom->rmass;
	double *vfrac = atom->vfrac;
	double *e = atom->e;
	double pInitial, d_iso, pFinal, p_rate, plastic_strain_increment;
	int i, itype;
	int nlocal = atom->nlocal;
	double dt = update->dt;
	double M_eff, p_wave_speed, mass_specific_energy, vol_specific_energy, rho;
	Matrix3d sigma_rate, eye, sigmaInitial, sigmaFinal, T, T_damaged, Jaumann_rate, sigma_rate_check;
	Matrix3d d_dev, sigmaInitial_dev, sigmaFinal_dev, sigma_dev_rate, strain;
	Vector3d x0i, xi, xp;

	eye.setIdentity();
	dtCFL = 1.0e22;
	pFinal = 0.0;

	for (i = 0; i < nlocal; i++) {
		particle_dt[i] = 0.0;

		itype = type[i];
		if (setflag[itype][itype] == 1) {
			if (mol[i] > 0) { // only do the following if particle has not failed -- mol < 0 means particle has failed

				/*
				 * initial stress state: given by the unrotateted Cauchy stress.
				 * Assemble Eigen 3d matrix from stored stress state
				 */
				sigmaInitial(0, 0) = tlsph_stress[i][0];
				sigmaInitial(0, 1) = tlsph_stress[i][1];
				sigmaInitial(0, 2) = tlsph_stress[i][2];
				sigmaInitial(1, 1) = tlsph_stress[i][3];
				sigmaInitial(1, 2) = tlsph_stress[i][4];
				sigmaInitial(2, 2) = tlsph_stress[i][5];
				sigmaInitial(1, 0) = sigmaInitial(0, 1);
				sigmaInitial(2, 0) = sigmaInitial(0, 2);
				sigmaInitial(2, 1) = sigmaInitial(1, 2);

				//cout << "this is sigma initial" << endl << sigmaInitial << endl;

				pInitial = sigmaInitial.trace() / 3.0; // isotropic part of initial stress
				sigmaInitial_dev = Deviator(sigmaInitial);
				d_iso = D[i].trace(); // volumetric part of stretch rate
				d_dev = Deviator(D[i]); // deviatoric part of stretch rate
				strain = 0.5 * (Fincr[i].transpose() * Fincr[i] - eye);
				mass_specific_energy = e[i] / rmass[i]; // energy per unit mass
				rho = rmass[i] / (detF[i] * vfrac[i]);
				vol_specific_energy = mass_specific_energy * rho; // energy per current volume

				/*
				 * pressure: compute pressure rate p_rate and final pressure pFinal
				 */

				ComputePressure(i, rho, mass_specific_energy, vol_specific_energy, pInitial, d_iso, pFinal, p_rate);

				/*
				 * material strength
				 */

				//cout << "this is the strain deviator rate" << endl << d_dev << endl;
				ComputeStressDeviator(i, sigmaInitial_dev, d_dev, sigmaFinal_dev, sigma_dev_rate, plastic_strain_increment);
				//cout << "this is the stress deviator rate" << endl << sigma_dev_rate << endl;

				// keep a rolling average of the plastic strain rate over the last 100 or so timesteps
				eff_plastic_strain[i] += plastic_strain_increment;

				// compute a characteristic time over which to average the plastic strain
				double tav = 1000 * radius[i] / (Lookup[SIGNAL_VELOCITY][itype]);
				eff_plastic_strain_rate[i] -= eff_plastic_strain_rate[i] / tav;
				eff_plastic_strain_rate[i] += (plastic_strain_increment / dt) / tav;
				eff_plastic_strain_rate[i] = MAX(0.0, eff_plastic_strain_rate[i]);

				/*
				 *  assemble total stress from pressure and deviatoric stress
				 */
				sigmaFinal = pFinal * eye + sigmaFinal_dev; // this is the stress that is kept

				if (JAUMANN) {
					/*
					 * sigma is already the co-rotated Cauchy stress.
					 * The stress rate, however, needs to be made objective.
					 */

					if (dt > 1.0e-16) {
						sigma_rate = (1.0 / dt) * (sigmaFinal - sigmaInitial);
					} else {
						sigma_rate.setZero();
					}

					Jaumann_rate = sigma_rate + W[i] * sigmaInitial + sigmaInitial * W[i].transpose();
					sigmaFinal = sigmaInitial + dt * Jaumann_rate;
					T = sigmaFinal;
				} else {
					/*
					 * sigma is the unrotated stress.
					 * need to do forward rotation of the unrotated stress sigma to the current configuration
					 */
					T = R[i] * sigmaFinal * R[i].transpose();
				}

				/*
				 * store unrotated stress in atom vector
				 * symmetry is exploited
				 */
				tlsph_stress[i][0] = sigmaFinal(0, 0);
				tlsph_stress[i][1] = sigmaFinal(0, 1);
				tlsph_stress[i][2] = sigmaFinal(0, 2);
				tlsph_stress[i][3] = sigmaFinal(1, 1);
				tlsph_stress[i][4] = sigmaFinal(1, 2);
				tlsph_stress[i][5] = sigmaFinal(2, 2);

				/*
				 *  Damage due to failure criteria.
				 */

				if (failureModel[itype].integration_point_wise) {
					ComputeDamage(i, strain, T, T_damaged);
					//T = T_damaged; Do not do this, it is undefined as of now
				}

				// store rotated, "true" Cauchy stress
				CauchyStress[i] = T;

				/*
				 * We have the corotational Cauchy stress.
				 * Convert to PK1. Note that reference configuration used for computing the forces is linked via
				 * the incremental deformation gradient, not the full deformation gradient.
				 */
				PK1[i] = detF[i] * T * FincrInv[i].transpose();

				/*
				 * pre-multiply stress tensor with shape matrix to save computation in force loop
				 */
				PK1[i] = PK1[i] * K[i];

				/*
				 * compute stable time step according to Pronto 2d
				 */

				Matrix3d deltaSigma;
				deltaSigma = sigmaFinal - sigmaInitial;
				p_rate = deltaSigma.trace() / (3.0 * dt + 1.0e-16);
				sigma_dev_rate = Deviator(deltaSigma) / (dt + 1.0e-16);

				double K_eff, mu_eff;
				effective_longitudinal_modulus(itype, dt, d_iso, p_rate, d_dev, sigma_dev_rate, damage[i], K_eff, mu_eff, M_eff);
				p_wave_speed = sqrt(M_eff / rho);

				if (mol[i] < 0) {
					error->one(FLERR, "this should not happen");
				}

				particle_dt[i] = 2.0 * radius[i] / p_wave_speed;
				dtCFL = MIN(dtCFL, particle_dt[i]);

			} else { // end if mol > 0
				PK1[i].setZero();
				K[i].setIdentity();
				CauchyStress[i].setZero();
				sigma_rate.setZero();
				tlsph_stress[i][0] = 0.0;
				tlsph_stress[i][1] = 0.0;
				tlsph_stress[i][2] = 0.0;
				tlsph_stress[i][3] = 0.0;
				tlsph_stress[i][4] = 0.0;
				tlsph_stress[i][5] = 0.0;
			} // end  if mol > 0
		} // end setflag
	} // end for
}

/* ----------------------------------------------------------------------
 allocate all arrays
 ------------------------------------------------------------------------- */

void PairTlsph::allocate() {
	allocated = 1;
	int n = atom->ntypes;

	memory->create(setflag, n + 1, n + 1, "pair:setflag");
	for (int i = 1; i <= n; i++)
		for (int j = i; j <= n; j++)
			setflag[i][j] = 0;

	memory->create(strengthModel, n + 1, "pair:strengthmodel");
	memory->create(eos, n + 1, "pair:eosmodel");
	failureModel = new failure_types[n + 1];
	memory->create(Lookup, MAX_KEY_VALUE, n + 1, "pair:LookupTable");

	memory->create(cutsq, n + 1, n + 1, "pair:cutsq"); // always needs to be allocated, even with granular neighborlist

	onerad_dynamic = new double[n + 1];
	onerad_frozen = new double[n + 1];
	maxrad_dynamic = new double[n + 1];
	maxrad_frozen = new double[n + 1];

}

/* ----------------------------------------------------------------------
 global settings
 ------------------------------------------------------------------------- */

void PairTlsph::settings(int narg, char **arg) {

	if (comm->me == 0) {
		printf(
				"\n>>========>>========>>========>>========>>========>>========>>========>>========>>========>>========>>========>>========\n");
		printf("TLSPH settings\n");
	}

	/*
	 * default value for update_threshold for updates of reference configuration:
	 * The maximum relative displacement which is tracked by the construction of LAMMPS' neighborlists
	 * is the folowing.
	 */

	cut_comm = MAX(neighbor->cutneighmax, comm->cutghostuser); // cutoff radius within which ghost atoms are communicated.
	update_threshold = cut_comm;
	update_method = UPDATE_NONE;

	int iarg = 0;

	while (true) {

		if (iarg >= narg) {
			break;
		}

		if (strcmp(arg[iarg], "*UPDATE_CONSTANT") == 0) {
			iarg++;
			if (iarg == narg) {
				error->all(FLERR, "expected number following *UPDATE_CONSTANT keyword");
			}

			update_method = UPDATE_CONSTANT_THRESHOLD;
			update_threshold = force->numeric(FLERR, arg[iarg]);

		} else if (strcmp(arg[iarg], "*UPDATE_PAIRWISE") == 0) {
			iarg++;
			if (iarg == narg) {
				error->all(FLERR, "expected number following *UPDATE_PAIRWISE keyword");
			}

			update_method = UPDATE_PAIRWISE_RATIO;
			update_threshold = force->numeric(FLERR, arg[iarg]);

		} else {
			char msg[128];
			sprintf(msg, "Illegal keyword for smd/integrate_tlsph: %s\n", arg[iarg]);
			error->all(FLERR, msg);
		}

		iarg++;
	}

	if ((update_threshold > cut_comm) && (update_method == UPDATE_CONSTANT_THRESHOLD)) {
		if (comm->me == 0) {
			printf("\n                ***** WARNING ***\n");
			printf("requested reference configuration update threshold is %g length units\n", update_threshold);
			printf("This value exceeds the maximum value %g beyond which TLSPH displacements can be tracked at current settings.\n",
					cut_comm);
			printf("Expect loss of neighbors!\n");
		}
	}

	if (comm->me == 0) {

		if (update_method == UPDATE_CONSTANT_THRESHOLD) {
			printf("... will update reference configuration if magnitude of relative displacement exceeds %g length units\n",
					update_threshold);
		} else if (update_method == UPDATE_PAIRWISE_RATIO) {
			printf("... will update reference configuration if ratio pairwise distance / smoothing length  exceeds %g\n",
					update_threshold);
		} else if (update_method == UPDATE_NONE) {
			printf("... will never update reference configuration");
		}
		printf(
				">>========>>========>>========>>========>>========>>========>>========>>========>>========>>========>>========>>========\n");

	}

}

/* ----------------------------------------------------------------------
 set coeffs for one or more type pairs
 ------------------------------------------------------------------------- */

void PairTlsph::coeff(int narg, char **arg) {
	int ioffset, iarg, iNextKwd, itype;
	char str[128];
	std::string s, t;

	if (narg < 3) {
		sprintf(str, "number of arguments for pair tlsph is too small!");
		error->all(FLERR, str);
	}
	if (!allocated)
		allocate();

	/*
	 * check that TLSPH parameters are given only in i,i form
	 */
	if (force->inumeric(FLERR, arg[0]) != force->inumeric(FLERR, arg[1])) {
		sprintf(str, "TLSPH coefficients can only be specified between particles of same type!");
		error->all(FLERR, str);
	}
	itype = force->inumeric(FLERR, arg[0]);

// set all eos, strength and failure models to inactive by default
	eos[itype] = EOS_NONE;
	strengthModel[itype] = STRENGTH_NONE;

	if (comm->me == 0) {
		printf(
				"\n>>========>>========>>========>>========>>========>>========>>========>>========>>========>>========>>========>>========\n");
		printf("SMD / TLSPH PROPERTIES OF PARTICLE TYPE %d:\n", itype);
	}

	/*
	 * read parameters which are common -- regardless of material / eos model
	 */

	ioffset = 2;
	if (strcmp(arg[ioffset], "*COMMON") != 0) {
		sprintf(str, "common keyword missing!");
		error->all(FLERR, str);
	}

	t = string("*");
	iNextKwd = -1;
	for (iarg = ioffset + 1; iarg < narg; iarg++) {
		s = string(arg[iarg]);
		if (s.compare(0, t.length(), t) == 0) {
			iNextKwd = iarg;
			break;
		}
	}

//printf("keyword following *COMMON is %s\n", arg[iNextKwd]);

	if (iNextKwd < 0) {
		sprintf(str, "no *KEYWORD terminates *COMMON");
		error->all(FLERR, str);
	}

	if (iNextKwd - ioffset != 7 + 1) {
		sprintf(str, "expected 7 arguments following *COMMON but got %d\n", iNextKwd - ioffset - 1);
		error->all(FLERR, str);
	}

	Lookup[REFERENCE_DENSITY][itype] = force->numeric(FLERR, arg[ioffset + 1]);
	Lookup[YOUNGS_MODULUS][itype] = force->numeric(FLERR, arg[ioffset + 2]);
	Lookup[POISSON_RATIO][itype] = force->numeric(FLERR, arg[ioffset + 3]);
	Lookup[VISCOSITY_Q1][itype] = force->numeric(FLERR, arg[ioffset + 4]);
	Lookup[VISCOSITY_Q2][itype] = force->numeric(FLERR, arg[ioffset + 5]);
	Lookup[HOURGLASS_CONTROL_AMPLITUDE][itype] = force->numeric(FLERR, arg[ioffset + 6]);
	Lookup[HEAT_CAPACITY][itype] = force->numeric(FLERR, arg[ioffset + 7]);

	Lookup[LAME_LAMBDA][itype] = Lookup[YOUNGS_MODULUS][itype] * Lookup[POISSON_RATIO][itype]
			/ ((1.0 + Lookup[POISSON_RATIO][itype] * (1.0 - 2.0 * Lookup[POISSON_RATIO][itype])));
	Lookup[SHEAR_MODULUS][itype] = Lookup[YOUNGS_MODULUS][itype] / (2.0 * (1.0 + Lookup[POISSON_RATIO][itype]));
	Lookup[M_MODULUS][itype] = Lookup[LAME_LAMBDA][itype] + 2.0 * Lookup[SHEAR_MODULUS][itype];
	Lookup[SIGNAL_VELOCITY][itype] = sqrt(
			(Lookup[LAME_LAMBDA][itype] + 2.0 * Lookup[SHEAR_MODULUS][itype]) / Lookup[REFERENCE_DENSITY][itype]);
	Lookup[BULK_MODULUS][itype] = Lookup[LAME_LAMBDA][itype] + 2.0 * Lookup[SHEAR_MODULUS][itype] / 3.0;

	if (comm->me == 0) {
		printf("\n material unspecific properties for SMD/TLSPH definition of particle type %d:\n", itype);
		printf("%60s : %g\n", "reference density", Lookup[REFERENCE_DENSITY][itype]);
		printf("%60s : %g\n", "Young's modulus", Lookup[YOUNGS_MODULUS][itype]);
		printf("%60s : %g\n", "Poisson ratio", Lookup[POISSON_RATIO][itype]);
		printf("%60s : %g\n", "linear viscosity coefficient", Lookup[VISCOSITY_Q1][itype]);
		printf("%60s : %g\n", "quadratic viscosity coefficient", Lookup[VISCOSITY_Q2][itype]);
		printf("%60s : %g\n", "hourglass control coefficient", Lookup[HOURGLASS_CONTROL_AMPLITUDE][itype]);
		printf("%60s : %g\n", "heat capacity [energy / (mass * temperature)]", Lookup[HEAT_CAPACITY][itype]);
		printf("%60s : %g\n", "Lame constant lambda", Lookup[LAME_LAMBDA][itype]);
		printf("%60s : %g\n", "shear modulus", Lookup[SHEAR_MODULUS][itype]);
		printf("%60s : %g\n", "bulk modulus", Lookup[BULK_MODULUS][itype]);
		printf("%60s : %g\n", "signal velocity", Lookup[SIGNAL_VELOCITY][itype]);

	}

	/*
	 * read following material cards
	 */

//printf("next kwd is %s\n", arg[iNextKwd]);
	eos[itype] = EOS_NONE;
	strengthModel[itype] = STRENGTH_NONE;

	while (true) {
		if (strcmp(arg[iNextKwd], "*END") == 0) {
			if (comm->me == 0) {
				printf("found *END keyword");
				printf(
						"\n>>========>>========>>========>>========>>========>>========>>========>>========>>========>>========>>========>>========\n\n");
			}
			break;
		}

		/*
		 * Linear Elasticity model based on deformation gradient
		 */
		ioffset = iNextKwd;
		if (strcmp(arg[ioffset], "*LINEAR_DEFGRAD") == 0) {
			strengthModel[itype] = LINEAR_DEFGRAD;

			if (comm->me == 0) {
				printf("reading *LINEAR_DEFGRAD\n");
			}

			t = string("*");
			iNextKwd = -1;
			for (iarg = ioffset + 1; iarg < narg; iarg++) {
				s = string(arg[iarg]);
				if (s.compare(0, t.length(), t) == 0) {
					iNextKwd = iarg;
					break;
				}
			}

			if (iNextKwd < 0) {
				sprintf(str, "no *KEYWORD terminates *LINEAR_DEFGRAD");
				error->all(FLERR, str);
			}

			if (iNextKwd - ioffset != 1) {
				sprintf(str, "expected 0 arguments following *LINEAR_DEFGRAD but got %d\n", iNextKwd - ioffset - 1);
				error->all(FLERR, str);
			}

			if (comm->me == 0) {
				printf("\n%60s\n", "Linear Elasticity model based on deformation gradient");
			}
		} else if (strcmp(arg[ioffset], "*STRENGTH_LINEAR") == 0) {

			/*
			 * Linear Elasticity strength only model based on strain rate
			 */

			strengthModel[itype] = STRENGTH_LINEAR;
			if (comm->me == 0) {
				printf("reading *STRENGTH_LINEAR\n");
			}

			t = string("*");
			iNextKwd = -1;
			for (iarg = ioffset + 1; iarg < narg; iarg++) {
				s = string(arg[iarg]);
				if (s.compare(0, t.length(), t) == 0) {
					iNextKwd = iarg;
					break;
				}
			}

			if (iNextKwd < 0) {
				sprintf(str, "no *KEYWORD terminates *STRENGTH_LINEAR");
				error->all(FLERR, str);
			}

			if (iNextKwd - ioffset != 1) {
				sprintf(str, "expected 0 arguments following *STRENGTH_LINEAR but got %d\n", iNextKwd - ioffset - 1);
				error->all(FLERR, str);
			}

			if (comm->me == 0) {
				printf("%60s\n", "Linear Elasticity strength based on strain rate");
			}
		} // end Linear Elasticity strength only model based on strain rate

		else if (strcmp(arg[ioffset], "*STRENGTH_LINEAR_PLASTIC") == 0) {

			/*
			 * Linear Elastic / perfectly plastic strength only model based on strain rate
			 */

			strengthModel[itype] = STRENGTH_LINEAR_PLASTIC;
			if (comm->me == 0) {
				printf("reading *STRENGTH_LINEAR_PLASTIC\n");
			}

			t = string("*");
			iNextKwd = -1;
			for (iarg = ioffset + 1; iarg < narg; iarg++) {
				s = string(arg[iarg]);
				if (s.compare(0, t.length(), t) == 0) {
					iNextKwd = iarg;
					break;
				}
			}

			if (iNextKwd < 0) {
				sprintf(str, "no *KEYWORD terminates *STRENGTH_LINEAR_PLASTIC");
				error->all(FLERR, str);
			}

			if (iNextKwd - ioffset != 2 + 1) {
				sprintf(str, "expected 2 arguments following *STRENGTH_LINEAR_PLASTIC but got %d\n", iNextKwd - ioffset - 1);
				error->all(FLERR, str);
			}

			Lookup[YIELD_STRESS][itype] = force->numeric(FLERR, arg[ioffset + 1]);
			Lookup[HARDENING_PARAMETER][itype] = force->numeric(FLERR, arg[ioffset + 2]);

			if (comm->me == 0) {
				printf("%60s\n", "Linear elastic / perfectly plastic strength based on strain rate");
				printf("%60s : %g\n", "Young's modulus", Lookup[YOUNGS_MODULUS][itype]);
				printf("%60s : %g\n", "Poisson ratio", Lookup[POISSON_RATIO][itype]);
				printf("%60s : %g\n", "shear modulus", Lookup[SHEAR_MODULUS][itype]);
				printf("%60s : %g\n", "constant yield stress", Lookup[YIELD_STRESS][itype]);
				printf("%60s : %g\n", "constant hardening parameter", Lookup[HARDENING_PARAMETER][itype]);
			}
		} // end Linear Elastic / perfectly plastic strength only model based on strain rate

		else if (strcmp(arg[ioffset], "*JOHNSON_COOK") == 0) {

			/*
			 * JOHNSON - COOK
			 */

			strengthModel[itype] = STRENGTH_JOHNSON_COOK;
			if (comm->me == 0) {
				printf("reading *JOHNSON_COOK\n");
			}

			t = string("*");
			iNextKwd = -1;
			for (iarg = ioffset + 1; iarg < narg; iarg++) {
				s = string(arg[iarg]);
				if (s.compare(0, t.length(), t) == 0) {
					iNextKwd = iarg;
					break;
				}
			}

			if (iNextKwd < 0) {
				sprintf(str, "no *KEYWORD terminates *JOHNSON_COOK");
				error->all(FLERR, str);
			}

			if (iNextKwd - ioffset != 8 + 1) {
				sprintf(str, "expected 8 arguments following *JOHNSON_COOK but got %d\n", iNextKwd - ioffset - 1);
				error->all(FLERR, str);
			}

			Lookup[JC_A][itype] = force->numeric(FLERR, arg[ioffset + 1]);
			Lookup[JC_B][itype] = force->numeric(FLERR, arg[ioffset + 2]);
			Lookup[JC_a][itype] = force->numeric(FLERR, arg[ioffset + 3]);
			Lookup[JC_C][itype] = force->numeric(FLERR, arg[ioffset + 4]);
			Lookup[JC_epdot0][itype] = force->numeric(FLERR, arg[ioffset + 5]);
			Lookup[JC_T0][itype] = force->numeric(FLERR, arg[ioffset + 6]);
			Lookup[JC_Tmelt][itype] = force->numeric(FLERR, arg[ioffset + 7]);
			Lookup[JC_M][itype] = force->numeric(FLERR, arg[ioffset + 8]);

			if (comm->me == 0) {
				printf("%60s\n", "Johnson Cook material strength model");
				printf("%60s : %g\n", "A: initial yield stress", Lookup[JC_A][itype]);
				printf("%60s : %g\n", "B : proportionality factor for plastic strain dependency", Lookup[JC_B][itype]);
				printf("%60s : %g\n", "a : exponent for plastic strain dependency", Lookup[JC_a][itype]);
				printf("%60s : %g\n", "C : proportionality factor for logarithmic plastic strain rate dependency",
						Lookup[JC_C][itype]);
				printf("%60s : %g\n", "epdot0 : dimensionality factor for plastic strain rate dependency",
						Lookup[JC_epdot0][itype]);
				printf("%60s : %g\n", "T0 : reference (room) temperature", Lookup[JC_T0][itype]);
				printf("%60s : %g\n", "Tmelt : melting temperature", Lookup[JC_Tmelt][itype]);
				printf("%60s : %g\n", "M : exponent for temperature dependency", Lookup[JC_M][itype]);
			}

		} else if (strcmp(arg[ioffset], "*EOS_NONE") == 0) {

			/*
			 * no eos
			 */

			eos[itype] = EOS_NONE;
			if (comm->me == 0) {
				printf("reading *EOS_NONE\n");
			}

			t = string("*");
			iNextKwd = -1;
			for (iarg = ioffset + 1; iarg < narg; iarg++) {
				s = string(arg[iarg]);
				if (s.compare(0, t.length(), t) == 0) {
					iNextKwd = iarg;
					break;
				}
			}

			if (iNextKwd < 0) {
				sprintf(str, "no *KEYWORD terminates *EOS_NONE");
				error->all(FLERR, str);
			}

			if (iNextKwd - ioffset != 1) {
				sprintf(str, "expected 0 arguments following *EOS_NONE but got %d\n", iNextKwd - ioffset - 1);
				error->all(FLERR, str);
			}

			if (comm->me == 0) {
				printf("\n%60s\n", "no EOS selected");
			}

		} else if (strcmp(arg[ioffset], "*EOS_LINEAR") == 0) {

			/*
			 * linear eos
			 */

			eos[itype] = EOS_LINEAR;
			if (comm->me == 0) {
				printf("reading *EOS_LINEAR\n");
			}

			t = string("*");
			iNextKwd = -1;
			for (iarg = ioffset + 1; iarg < narg; iarg++) {
				s = string(arg[iarg]);
				if (s.compare(0, t.length(), t) == 0) {
					iNextKwd = iarg;
					break;
				}
			}

			if (iNextKwd < 0) {
				sprintf(str, "no *KEYWORD terminates *EOS_LINEAR");
				error->all(FLERR, str);
			}

			if (iNextKwd - ioffset != 1) {
				sprintf(str, "expected 0 arguments following *EOS_LINEAR but got %d\n", iNextKwd - ioffset - 1);
				error->all(FLERR, str);
			}

			if (comm->me == 0) {
				printf("\n%60s\n", "linear EOS based on strain rate");
				printf("%60s : %g\n", "bulk modulus", Lookup[BULK_MODULUS][itype]);
			}
		} // end linear eos
		else if (strcmp(arg[ioffset], "*EOS_SHOCK") == 0) {

			/*
			 * shock eos
			 */

			eos[itype] = EOS_SHOCK;
			if (comm->me == 0) {
				printf("reading *EOS_SHOCK\n");
			}

			t = string("*");
			iNextKwd = -1;
			for (iarg = ioffset + 1; iarg < narg; iarg++) {
				s = string(arg[iarg]);
				if (s.compare(0, t.length(), t) == 0) {
					iNextKwd = iarg;
					break;
				}
			}

			if (iNextKwd < 0) {
				sprintf(str, "no *KEYWORD terminates *EOS_SHOCK");
				error->all(FLERR, str);
			}

			if (iNextKwd - ioffset != 3 + 1) {
				sprintf(str, "expected 3 arguments (c0, S, Gamma) following *EOS_SHOCK but got %d\n", iNextKwd - ioffset - 1);
				error->all(FLERR, str);
			}

			Lookup[EOS_SHOCK_C0][itype] = force->numeric(FLERR, arg[ioffset + 1]);
			Lookup[EOS_SHOCK_S][itype] = force->numeric(FLERR, arg[ioffset + 2]);
			Lookup[EOS_SHOCK_GAMMA][itype] = force->numeric(FLERR, arg[ioffset + 3]);
			if (comm->me == 0) {
				printf("\n%60s\n", "shock EOS based on strain rate");
				printf("%60s : %g\n", "reference speed of sound", Lookup[EOS_SHOCK_C0][itype]);
				printf("%60s : %g\n", "Hugoniot parameter S", Lookup[EOS_SHOCK_S][itype]);
				printf("%60s : %g\n", "Grueneisen Gamma", Lookup[EOS_SHOCK_GAMMA][itype]);
			}
		} // end shock eos

		else if (strcmp(arg[ioffset], "*EOS_POLYNOMIAL") == 0) {
			/*
			 * polynomial eos
			 */

			eos[itype] = EOS_POLYNOMIAL;
			if (comm->me == 0) {
				printf("reading *EOS_POLYNOMIAL\n");
			}

			t = string("*");
			iNextKwd = -1;
			for (iarg = ioffset + 1; iarg < narg; iarg++) {
				s = string(arg[iarg]);
				if (s.compare(0, t.length(), t) == 0) {
					iNextKwd = iarg;
					break;
				}
			}

			if (iNextKwd < 0) {
				sprintf(str, "no *KEYWORD terminates *EOS_POLYNOMIAL");
				error->all(FLERR, str);
			}

			if (iNextKwd - ioffset != 7 + 1) {
				sprintf(str, "expected 7 arguments following *EOS_POLYNOMIAL but got %d\n", iNextKwd - ioffset - 1);
				error->all(FLERR, str);
			}

			Lookup[EOS_POLYNOMIAL_C0][itype] = force->numeric(FLERR, arg[ioffset + 1]);
			Lookup[EOS_POLYNOMIAL_C1][itype] = force->numeric(FLERR, arg[ioffset + 2]);
			Lookup[EOS_POLYNOMIAL_C2][itype] = force->numeric(FLERR, arg[ioffset + 3]);
			Lookup[EOS_POLYNOMIAL_C3][itype] = force->numeric(FLERR, arg[ioffset + 4]);
			Lookup[EOS_POLYNOMIAL_C4][itype] = force->numeric(FLERR, arg[ioffset + 5]);
			Lookup[EOS_POLYNOMIAL_C5][itype] = force->numeric(FLERR, arg[ioffset + 6]);
			Lookup[EOS_POLYNOMIAL_C6][itype] = force->numeric(FLERR, arg[ioffset + 7]);
			if (comm->me == 0) {
				printf("\n%60s\n", "polynomial EOS based on strain rate");
				printf("%60s : %g\n", "parameter c0", Lookup[EOS_POLYNOMIAL_C0][itype]);
				printf("%60s : %g\n", "parameter c1", Lookup[EOS_POLYNOMIAL_C1][itype]);
				printf("%60s : %g\n", "parameter c2", Lookup[EOS_POLYNOMIAL_C2][itype]);
				printf("%60s : %g\n", "parameter c3", Lookup[EOS_POLYNOMIAL_C3][itype]);
				printf("%60s : %g\n", "parameter c4", Lookup[EOS_POLYNOMIAL_C4][itype]);
				printf("%60s : %g\n", "parameter c5", Lookup[EOS_POLYNOMIAL_C5][itype]);
				printf("%60s : %g\n", "parameter c6", Lookup[EOS_POLYNOMIAL_C6][itype]);
			}
		} // end polynomial eos

		else if (strcmp(arg[ioffset], "*FAILURE_MAX_PLASTIC_STRAIN") == 0) {

			/*
			 * maximum plastic strain failure criterion
			 */

			if (comm->me == 0) {
				printf("reading *FAILURE_MAX_PLASTIC_SRTRAIN\n");
			}

			t = string("*");
			iNextKwd = -1;
			for (iarg = ioffset + 1; iarg < narg; iarg++) {
				s = string(arg[iarg]);
				if (s.compare(0, t.length(), t) == 0) {
					iNextKwd = iarg;
					break;
				}
			}

			if (iNextKwd < 0) {
				sprintf(str, "no *KEYWORD terminates *FAILURE_MAX_PLASTIC_STRAIN");
				error->all(FLERR, str);
			}

			if (iNextKwd - ioffset != 1 + 1) {
				sprintf(str, "expected 1 arguments following *FAILURE_MAX_PLASTIC_STRAIN but got %d\n", iNextKwd - ioffset - 1);
				error->all(FLERR, str);
			}

			failureModel[itype].failure_max_plastic_strain = true;
			failureModel[itype].integration_point_wise = true;
			Lookup[FAILURE_MAX_PLASTIC_STRAIN_THRESHOLD][itype] = force->numeric(FLERR, arg[ioffset + 1]);

			if (comm->me == 0) {
				printf("\n%60s\n", "maximum plastic strain failure criterion");
				printf("%60s : %g\n", "failure occurs when plastic strain reaches limit",
						Lookup[FAILURE_MAX_PLASTIC_STRAIN_THRESHOLD][itype]);
			}
		} // end maximum plastic strain failure criterion
		else if (strcmp(arg[ioffset], "*FAILURE_MAX_PAIRWISE_STRAIN") == 0) {

			/*
			 * failure criterion based on maximum strain between a pair of TLSPH particles.
			 */

			if (comm->me == 0) {
				printf("reading *FAILURE_MAX_PAIRWISE_STRAIN\n");
			}

			if (update_method != UPDATE_NONE) {
				error->all(FLERR, "cannot use *FAILURE_MAX_PAIRWISE_STRAIN with updated Total-Lagrangian formalism");
			}

			t = string("*");
			iNextKwd = -1;
			for (iarg = ioffset + 1; iarg < narg; iarg++) {
				s = string(arg[iarg]);
				if (s.compare(0, t.length(), t) == 0) {
					iNextKwd = iarg;
					break;
				}
			}

			if (iNextKwd < 0) {
				sprintf(str, "no *KEYWORD terminates *FAILURE_MAX_PAIRWISE_STRAIN");
				error->all(FLERR, str);
			}

			if (iNextKwd - ioffset != 1 + 1) {
				sprintf(str, "expected 1 arguments following *FAILURE_MAX_PAIRWISE_STRAIN but got %d\n", iNextKwd - ioffset - 1);
				error->all(FLERR, str);
			}

			failureModel[itype].failure_max_pairwise_strain = true;
			failureModel[itype].integration_point_wise = true;
			Lookup[FAILURE_MAX_PAIRWISE_STRAIN_THRESHOLD][itype] = force->numeric(FLERR, arg[ioffset + 1]);

			if (comm->me == 0) {
				printf("\n%60s\n", "maximum pairwise strain failure criterion");
				printf("%60s : %g\n", "failure occurs when pairwise strain reaches limit",
						Lookup[FAILURE_MAX_PAIRWISE_STRAIN_THRESHOLD][itype]);
			}
		} // end pair based maximum strain failure criterion
		else if (strcmp(arg[ioffset], "*FAILURE_MAX_PRINCIPAL_STRAIN") == 0) {
			error->all(FLERR, "this failure model is currently unsupported");

			/*
			 * maximum principal strain failure criterion
			 */
			if (comm->me == 0) {
				printf("reading *FAILURE_MAX_PRINCIPAL_STRAIN\n");
			}

			t = string("*");
			iNextKwd = -1;
			for (iarg = ioffset + 1; iarg < narg; iarg++) {
				s = string(arg[iarg]);
				if (s.compare(0, t.length(), t) == 0) {
					iNextKwd = iarg;
					break;
				}
			}

			if (iNextKwd < 0) {
				sprintf(str, "no *KEYWORD terminates *FAILURE_MAX_PRINCIPAL_STRAIN");
				error->all(FLERR, str);
			}

			if (iNextKwd - ioffset != 1 + 1) {
				sprintf(str, "expected 1 arguments following *FAILURE_MAX_PRINCIPAL_STRAIN but got %d\n", iNextKwd - ioffset - 1);
				error->all(FLERR, str);
			}

			failureModel[itype].failure_max_principal_strain = true;
			failureModel[itype].integration_point_wise = true;
			Lookup[FAILURE_MAX_PRINCIPAL_STRAIN_THRESHOLD][itype] = force->numeric(FLERR, arg[ioffset + 1]);

			if (comm->me == 0) {
				printf("\n%60s\n", "maximum principal strain failure criterion");
				printf("%60s : %g\n", "failure occurs when principal strain reaches limit",
						Lookup[FAILURE_MAX_PRINCIPAL_STRAIN_THRESHOLD][itype]);
			}
		} // end maximum principal strain failure criterion
		else if (strcmp(arg[ioffset], "*FAILURE_JOHNSON_COOK") == 0) {
			error->all(FLERR, "this failure model is currently unsupported");
			if (comm->me == 0) {
				printf("reading *FAILURE_JOHNSON_COOK\n");
			}

			t = string("*");
			iNextKwd = -1;
			for (iarg = ioffset + 1; iarg < narg; iarg++) {
				s = string(arg[iarg]);
				if (s.compare(0, t.length(), t) == 0) {
					iNextKwd = iarg;
					break;
				}
			}

			if (iNextKwd < 0) {
				sprintf(str, "no *KEYWORD terminates *FAILURE_JOHNSON_COOK");
				error->all(FLERR, str);
			}

			if (iNextKwd - ioffset != 5 + 1) {
				sprintf(str, "expected 5 arguments following *FAILURE_JOHNSON_COOK but got %d\n", iNextKwd - ioffset - 1);
				error->all(FLERR, str);
			}

			failureModel[itype].failure_johnson_cook = true;
			failureModel[itype].integration_point_wise = true;

			Lookup[FAILURE_JC_D1][itype] = force->numeric(FLERR, arg[ioffset + 1]);
			Lookup[FAILURE_JC_D2][itype] = force->numeric(FLERR, arg[ioffset + 2]);
			Lookup[FAILURE_JC_D3][itype] = force->numeric(FLERR, arg[ioffset + 3]);
			Lookup[FAILURE_JC_D4][itype] = force->numeric(FLERR, arg[ioffset + 4]);
			Lookup[FAILURE_JC_EPDOT0][itype] = force->numeric(FLERR, arg[ioffset + 5]);

			if (comm->me == 0) {
				printf("\n%60s\n", "Johnson-Cook failure criterion");
				printf("%60s : %g\n", "parameter d1", Lookup[FAILURE_JC_D1][itype]);
				printf("%60s : %g\n", "parameter d2", Lookup[FAILURE_JC_D2][itype]);
				printf("%60s : %g\n", "parameter d3", Lookup[FAILURE_JC_D3][itype]);
				printf("%60s : %g\n", "parameter d4", Lookup[FAILURE_JC_D4][itype]);
				printf("%60s : %g\n", "reference plastic strain rate", Lookup[FAILURE_JC_EPDOT0][itype]);
			}

		} else if (strcmp(arg[ioffset], "*FAILURE_MAX_PRINCIPAL_STRESS") == 0) {
			error->all(FLERR, "this failure model is currently unsupported");

			/*
			 * maximum principal stress failure criterion
			 */

			if (comm->me == 0) {
				printf("reading *FAILURE_MAX_PRINCIPAL_STRESS\n");
			}

			t = string("*");
			iNextKwd = -1;
			for (iarg = ioffset + 1; iarg < narg; iarg++) {
				s = string(arg[iarg]);
				if (s.compare(0, t.length(), t) == 0) {
					iNextKwd = iarg;
					break;
				}
			}

			if (iNextKwd < 0) {
				sprintf(str, "no *KEYWORD terminates *FAILURE_MAX_PRINCIPAL_STRESS");
				error->all(FLERR, str);
			}

			if (iNextKwd - ioffset != 1 + 1) {
				sprintf(str, "expected 1 arguments following *FAILURE_MAX_PRINCIPAL_STRESS but got %d\n", iNextKwd - ioffset - 1);
				error->all(FLERR, str);
			}

			failureModel[itype].failure_max_principal_stress = true;
			failureModel[itype].integration_point_wise = true;
			Lookup[FAILURE_MAX_PRINCIPAL_STRESS_THRESHOLD][itype] = force->numeric(FLERR, arg[ioffset + 1]);

			if (comm->me == 0) {
				printf("\n%60s\n", "maximum principal stress failure criterion");
				printf("%60s : %g\n", "failure occurs when principal stress reaches limit",
						Lookup[FAILURE_MAX_PRINCIPAL_STRESS_THRESHOLD][itype]);
			}
		} // end maximum principal stress failure criterion

		else if (strcmp(arg[ioffset], "*FAILURE_ENERGY_RELEASE_RATE") == 0) {
			if (comm->me == 0) {
				printf("reading *FAILURE_ENERGY_RELEASE_RATE\n");
			}

			t = string("*");
			iNextKwd = -1;
			for (iarg = ioffset + 1; iarg < narg; iarg++) {
				s = string(arg[iarg]);
				if (s.compare(0, t.length(), t) == 0) {
					iNextKwd = iarg;
					break;
				}
			}

			if (iNextKwd < 0) {
				sprintf(str, "no *KEYWORD terminates *FAILURE_ENERGY_RELEASE_RATE");
				error->all(FLERR, str);
			}

			if (iNextKwd - ioffset != 1 + 1) {
				sprintf(str, "expected 1 arguments following *FAILURE_ENERGY_RELEASE_RATE but got %d\n", iNextKwd - ioffset - 1);
				error->all(FLERR, str);
			}

			failureModel[itype].failure_energy_release_rate = true;
			Lookup[CRITICAL_ENERGY_RELEASE_RATE][itype] = force->numeric(FLERR, arg[ioffset + 1]);

			if (comm->me == 0) {
				printf("\n%60s\n", "critical energy release rate failure criterion");
				printf("%60s : %g\n", "failure occurs when energy release rate reaches limit",
						Lookup[CRITICAL_ENERGY_RELEASE_RATE][itype]);
			}
		} // end energy release rate failure criterion

		else {
			sprintf(str, "unknown *KEYWORD: %s", arg[ioffset]);
			error->all(FLERR, str);
		}

	}

	setflag[itype][itype] = 1;

}

/* ----------------------------------------------------------------------
 init for one type pair i,j and corresponding j,i
 ------------------------------------------------------------------------- */

double PairTlsph::init_one(int i, int j) {

	if (!allocated)
		allocate();

	if (setflag[i][j] == 0)
		error->all(FLERR, "All pair coeffs are not set");

	if (force->newton == 1)
		error->all(FLERR, "Pair style tlsph requires newton off");

// cutoff = sum of max I,J radii for
// dynamic/dynamic & dynamic/frozen interactions, but not frozen/frozen

	double cutoff = maxrad_dynamic[i] + maxrad_dynamic[j];
	cutoff = MAX(cutoff, maxrad_frozen[i] + maxrad_dynamic[j]);
	cutoff = MAX(cutoff, maxrad_dynamic[i] + maxrad_frozen[j]);
//printf("cutoff for pair pair tlsph = %f\n", cutoff);
	return cutoff;
}

/* ----------------------------------------------------------------------
 init specific to this pair style
 ------------------------------------------------------------------------- */

void PairTlsph::init_style() {
	int i;

	if (force->newton_pair == 1) {
		error->all(FLERR, "Pair style tlsph requires newton pair off");
	}

// request a granular neighbor list
	int irequest = neighbor->request(this);
	neighbor->requests[irequest]->half = 0;
	neighbor->requests[irequest]->gran = 1;

// set maxrad_dynamic and maxrad_frozen for each type
// include future Fix pour particles as dynamic

	for (i = 1; i <= atom->ntypes; i++)
		onerad_dynamic[i] = onerad_frozen[i] = 0.0;

	double *radius = atom->radius;
	int *type = atom->type;
	int nlocal = atom->nlocal;

	for (i = 0; i < nlocal; i++)
		onerad_dynamic[type[i]] = MAX(onerad_dynamic[type[i]], radius[i]);

	MPI_Allreduce(&onerad_dynamic[1], &maxrad_dynamic[1], atom->ntypes, MPI_DOUBLE, MPI_MAX, world);
	MPI_Allreduce(&onerad_frozen[1], &maxrad_frozen[1], atom->ntypes, MPI_DOUBLE, MPI_MAX, world);

// if first init, create Fix needed for storing reference configuration neighbors

	int igroup = group->find("tlsph");
	if (igroup == -1)
		error->all(FLERR, "Pair style tlsph requires its particles to be part of a group named tlsph. This group does not exist.");

	if (fix_tlsph_reference_configuration == NULL) {
		char **fixarg = new char*[3];
		fixarg[0] = (char *) "SMD_TLSPH_NEIGHBORS";
		fixarg[1] = (char *) "tlsph";
		fixarg[2] = (char *) "SMD_TLSPH_NEIGHBORS";
		modify->add_fix(3, fixarg);
		delete[] fixarg;
		fix_tlsph_reference_configuration = (FixSMD_TLSPH_ReferenceConfiguration *) modify->fix[modify->nfix - 1];
		fix_tlsph_reference_configuration->pair = this;
	}

// find associated SMD_TLSPH_NEIGHBORS fix that must exist
// could have changed locations in fix list since created

	ifix_tlsph = -1;
	for (int i = 0; i < modify->nfix; i++)
		if (strcmp(modify->fix[i]->style, "SMD_TLSPH_NEIGHBORS") == 0)
			ifix_tlsph = i;
	if (ifix_tlsph == -1)
		error->all(FLERR, "Fix SMD_TLSPH_NEIGHBORS does not exist");

}

/* ----------------------------------------------------------------------
 neighbor callback to inform pair style of neighbor list to use
 optional granular history list
 ------------------------------------------------------------------------- */

void PairTlsph::init_list(int id, NeighList *ptr) {
	if (id == 0)
		list = ptr;
}

/* ----------------------------------------------------------------------
 memory usage of local atom-based arrays
 ------------------------------------------------------------------------- */

double PairTlsph::memory_usage() {

	return 118 * nmax * sizeof(double);
}

/* ----------------------------------------------------------------------
 extract method to provide access to this class' data structures
 ------------------------------------------------------------------------- */

void *PairTlsph::extract(const char *str, int &i) {
//printf("in PairTlsph::extract\n");
	if (strcmp(str, "smd/tlsph/Fincr_ptr") == 0) {
		return (void *) Fincr;
	} else if (strcmp(str, "smd/tlsph/detF_ptr") == 0) {
		return (void *) detF;
	} else if (strcmp(str, "smd/tlsph/PK1_ptr") == 0) {
		return (void *) PK1;
	} else if (strcmp(str, "smd/tlsph/smoothVel_ptr") == 0) {
		return (void *) smoothVelDifference;
	} else if (strcmp(str, "smd/tlsph/numNeighsRefConfig_ptr") == 0) {
		return (void *) numNeighsRefConfig;
	} else if (strcmp(str, "smd/tlsph/stressTensor_ptr") == 0) {
		return (void *) CauchyStress;
	} else if (strcmp(str, "smd/tlsph/updateFlag_ptr") == 0) {
		return (void *) &updateFlag;
	} else if (strcmp(str, "smd/tlsph/strain_rate_ptr") == 0) {
		return (void *) D;
	} else if (strcmp(str, "smd/tlsph/hMin_ptr") == 0) {
		return (void *) &hMin;
	} else if (strcmp(str, "smd/tlsph/dtCFL_ptr") == 0) {
		return (void *) &dtCFL;
	} else if (strcmp(str, "smd/tlsph/dtRelative_ptr") == 0) {
		return (void *) &dtRelative;
	} else if (strcmp(str, "smd/tlsph/hourglass_error_ptr") == 0) {
		return (void *) hourglass_error;
	} else if (strcmp(str, "smd/tlsph/particle_dt_ptr") == 0) {
		return (void *) particle_dt;
	} else if (strcmp(str, "smd/tlsph/rotation_ptr") == 0) {
		return (void *) R;
	}

	return NULL;
}

/* ---------------------------------------------------------------------- */

int PairTlsph::pack_forward_comm(int n, int *list, double *buf, int pbc_flag, int *pbc) {
	int i, j, m;
	tagint *mol = atom->molecule;
	double *damage = atom->damage;
	double *eff_plastic_strain = atom->eff_plastic_strain;
	double *eff_plastic_strain_rate = atom->eff_plastic_strain_rate;

//printf("in PairTlsph::pack_forward_comm\n");

	m = 0;
	for (i = 0; i < n; i++) {
		j = list[i];
		buf[m++] = PK1[j](0, 0); // PK1 is not symmetric
		buf[m++] = PK1[j](0, 1);
		buf[m++] = PK1[j](0, 2);
		buf[m++] = PK1[j](1, 0);
		buf[m++] = PK1[j](1, 1);
		buf[m++] = PK1[j](1, 2);
		buf[m++] = PK1[j](2, 0);
		buf[m++] = PK1[j](2, 1);
		buf[m++] = PK1[j](2, 2); // 9

		buf[m++] = Fincr[j](0, 0); // Fincr is not symmetric
		buf[m++] = Fincr[j](0, 1);
		buf[m++] = Fincr[j](0, 2);
		buf[m++] = Fincr[j](1, 0);
		buf[m++] = Fincr[j](1, 1);
		buf[m++] = Fincr[j](1, 2);
		buf[m++] = Fincr[j](2, 0);
		buf[m++] = Fincr[j](2, 1);
		buf[m++] = Fincr[j](2, 2); // 9 + 9 = 18

		buf[m++] = mol[j]; //19
		buf[m++] = damage[j]; //20
		buf[m++] = eff_plastic_strain[j]; //21
		buf[m++] = eff_plastic_strain_rate[j]; //22

	}
	return m;
}

/* ---------------------------------------------------------------------- */

void PairTlsph::unpack_forward_comm(int n, int first, double *buf) {
	int i, m, last;
	tagint *mol = atom->molecule;
	double *damage = atom->damage;
	double *eff_plastic_strain = atom->eff_plastic_strain;
	double *eff_plastic_strain_rate = atom->eff_plastic_strain_rate;

//printf("in PairTlsph::unpack_forward_comm\n");

	m = 0;
	last = first + n;
	for (i = first; i < last; i++) {

		PK1[i](0, 0) = buf[m++]; // PK1 is not symmetric
		PK1[i](0, 1) = buf[m++];
		PK1[i](0, 2) = buf[m++];
		PK1[i](1, 0) = buf[m++];
		PK1[i](1, 1) = buf[m++];
		PK1[i](1, 2) = buf[m++];
		PK1[i](2, 0) = buf[m++];
		PK1[i](2, 1) = buf[m++];
		PK1[i](2, 2) = buf[m++];

		Fincr[i](0, 0) = buf[m++];
		Fincr[i](0, 1) = buf[m++];
		Fincr[i](0, 2) = buf[m++];
		Fincr[i](1, 0) = buf[m++];
		Fincr[i](1, 1) = buf[m++];
		Fincr[i](1, 2) = buf[m++];
		Fincr[i](2, 0) = buf[m++];
		Fincr[i](2, 1) = buf[m++];
		Fincr[i](2, 2) = buf[m++];

		mol[i] = static_cast<int>(buf[m++]);
		damage[i] = buf[m++];
		eff_plastic_strain[i] = buf[m++]; //22
		eff_plastic_strain_rate[i] = buf[m++]; //23
	}
}

/* ----------------------------------------------------------------------
 compute effective P-wave speed
 determined by longitudinal modulus
 ------------------------------------------------------------------------- */

void PairTlsph::effective_longitudinal_modulus(const int itype, const double dt, const double d_iso, const double p_rate,
		const Matrix3d d_dev, const Matrix3d sigma_dev_rate, const double damage, double &K_eff, double &mu_eff, double &M_eff) {
	double M0; // initial longitudinal modulus
	double shear_rate_sq;

//	if (damage >= 0.5) {
//		M_eff = Lookup[M_MODULUS][itype];
//		K_eff = Lookup[BULK_MODULUS][itype];
//		mu_eff = Lookup[SHEAR_MODULUS][itype];
//		return;
//	}

	M0 = Lookup[M_MODULUS][itype];

	if (dt * d_iso > 1.0e-6) {
		K_eff = p_rate / d_iso;
		if (K_eff < 0.0) { // it is possible for K_eff to become negative due to strain softening
//			if (damage == 0.0) {
//				error->one(FLERR, "computed a negative effective bulk modulus but particle is not damaged.");
//			}
			K_eff = Lookup[BULK_MODULUS][itype];
		}
	} else {
		K_eff = Lookup[BULK_MODULUS][itype];
	}

	if (domain->dimension == 3) {
// Calculate 2 mu by looking at ratio shear stress / shear strain. Use numerical softening to avoid divide-by-zero.
		mu_eff = 0.5
				* (sigma_dev_rate(0, 1) / (d_dev(0, 1) + 1.0e-16) + sigma_dev_rate(0, 2) / (d_dev(0, 2) + 1.0e-16)
						+ sigma_dev_rate(1, 2) / (d_dev(1, 2) + 1.0e-16));

// Calculate magnitude of deviatoric strain rate. This is used for deciding if shear modulus should be computed from current rate or be taken as the initial value.
		shear_rate_sq = d_dev(0, 1) * d_dev(0, 1) + d_dev(0, 2) * d_dev(0, 2) + d_dev(1, 2) * d_dev(1, 2);
	} else {
		mu_eff = 0.5 * (sigma_dev_rate(0, 1) / (d_dev(0, 1) + 1.0e-16));
		shear_rate_sq = d_dev(0, 1) * d_dev(0, 1);
	}

	if (dt * dt * shear_rate_sq < 1.0e-8) {
		mu_eff = Lookup[SHEAR_MODULUS][itype];
	}

	if (mu_eff < Lookup[SHEAR_MODULUS][itype]) { // it is possible for mu_eff to become negative due to strain softening
//		if (damage == 0.0) {
//			printf("mu_eff = %f, tau=%f, gamma=%f\n", mu_eff, sigma_dev_rate(0, 1), d_dev(0, 1));
//			error->message(FLERR, "computed a negative effective shear modulus but particle is not damaged.");
//		}
		mu_eff = Lookup[SHEAR_MODULUS][itype];
	}

//mu_eff = Lookup[SHEAR_MODULUS][itype];

	if (K_eff < 0.0) {
		printf("K_eff = %f, p_rate=%f, vol_rate=%f\n", K_eff, p_rate, d_iso);
	}

	if (mu_eff < 0.0) {
		printf("mu_eff = %f, tau=%f, gamma=%f\n", mu_eff, sigma_dev_rate(0, 1), d_dev(0, 1));
		error->one(FLERR, "");
	}

	M_eff = (K_eff + 4.0 * mu_eff / 3.0); // effective dilational modulus, see Pronto 2d eqn 3.4.8

	if (M_eff < M0) { // do not allow effective dilatational modulus to decrease beyond its initial value
		M_eff = M0;
	}
}

/* ----------------------------------------------------------------------
 compute pressure. Called from AssembleStress().
 ------------------------------------------------------------------------- */
void PairTlsph::ComputePressure(const int i, const double rho, const double mass_specific_energy, const double vol_specific_energy,
	const double pInitial, const double d_iso, double &pFinal, double &p_rate) {
	int *type = atom->type;
	double dt = update->dt;

	int itype;

	itype = type[i];

	switch (eos[itype]) {
	case EOS_LINEAR:
		LinearEOS(Lookup[BULK_MODULUS][itype], pInitial, d_iso, dt, pFinal, p_rate);
		break;
	case EOS_NONE:
		pFinal = 0.0;
		p_rate = 0.0;
		break;
	case EOS_SHOCK:
//  rho,  rho0,  e,  e0,  c0,  S,  Gamma,  pInitial,  dt,  &pFinal,  &p_rate);
		ShockEOS(rho, Lookup[REFERENCE_DENSITY][itype], mass_specific_energy, 0.0, Lookup[EOS_SHOCK_C0][itype],
				Lookup[EOS_SHOCK_S][itype], Lookup[EOS_SHOCK_GAMMA][itype], pInitial, dt, pFinal, p_rate);
		break;
	case EOS_POLYNOMIAL:
		polynomialEOS(rho, Lookup[REFERENCE_DENSITY][itype], vol_specific_energy, Lookup[EOS_POLYNOMIAL_C0][itype],
				Lookup[EOS_POLYNOMIAL_C1][itype], Lookup[EOS_POLYNOMIAL_C2][itype], Lookup[EOS_POLYNOMIAL_C3][itype],
				Lookup[EOS_POLYNOMIAL_C4][itype], Lookup[EOS_POLYNOMIAL_C5][itype], Lookup[EOS_POLYNOMIAL_C6][itype], pInitial, dt,
				pFinal, p_rate);

		break;
	default:
		error->one(FLERR, "unknown EOS.");
		break;
	}
}

/* ----------------------------------------------------------------------
 Compute stress deviator. Called from AssembleStress().
 ------------------------------------------------------------------------- */
void PairTlsph::ComputeStressDeviator(const int i, const Matrix3d sigmaInitial_dev, const Matrix3d d_dev, Matrix3d &sigmaFinal_dev,
		Matrix3d &sigma_dev_rate, double &plastic_strain_increment) {
	double *eff_plastic_strain = atom->eff_plastic_strain;
	double *eff_plastic_strain_rate = atom->eff_plastic_strain_rate;
	int *type = atom->type;
	double *rmass = atom->rmass;
//double *vfrac = atom->vfrac;
	double *e = atom->e;
	double dt = update->dt;
	double yieldStress;
	int itype;

	double mass_specific_energy = e[i] / rmass[i]; // energy per unit mass
	plastic_strain_increment = 0.0;
	itype = type[i];

	switch (strengthModel[itype]) {
	case STRENGTH_LINEAR:

		sigma_dev_rate = 2.0 * Lookup[SHEAR_MODULUS][itype] * d_dev;
		sigmaFinal_dev = sigmaInitial_dev + dt * sigma_dev_rate;

		break;
	case LINEAR_DEFGRAD:
//LinearStrengthDefgrad(Lookup[LAME_LAMBDA][itype], Lookup[SHEAR_MODULUS][itype], Fincr[i], &sigmaFinal_dev);
//eff_plastic_strain[i] = 0.0;
//p_rate = pInitial - sigmaFinal_dev.trace() / 3.0;
//sigma_dev_rate = sigmaInitial_dev - Deviator(sigmaFinal_dev);
		error->one(FLERR, "LINEAR_DEFGRAD is only for debugging purposes and currently deactivated.");
		R[i].setIdentity();
		break;
	case STRENGTH_LINEAR_PLASTIC:

		yieldStress = Lookup[YIELD_STRESS][itype] + Lookup[HARDENING_PARAMETER][itype] * eff_plastic_strain[i];
		LinearPlasticStrength(Lookup[SHEAR_MODULUS][itype], yieldStress, sigmaInitial_dev, d_dev, dt, sigmaFinal_dev,
				sigma_dev_rate, plastic_strain_increment);
		break;
	case STRENGTH_JOHNSON_COOK:
		JohnsonCookStrength(Lookup[SHEAR_MODULUS][itype], Lookup[HEAT_CAPACITY][itype], mass_specific_energy, Lookup[JC_A][itype],
				Lookup[JC_B][itype], Lookup[JC_a][itype], Lookup[JC_C][itype], Lookup[JC_epdot0][itype], Lookup[JC_T0][itype],
				Lookup[JC_Tmelt][itype], Lookup[JC_M][itype], dt, eff_plastic_strain[i], eff_plastic_strain_rate[i],
				sigmaInitial_dev, d_dev, sigmaFinal_dev, sigma_dev_rate, plastic_strain_increment);
		break;
	case STRENGTH_NONE:
		sigmaFinal_dev.setZero();
		sigma_dev_rate.setZero();
		break;
	default:
		error->one(FLERR, "unknown strength model.");
		break;
	}

}

/* ----------------------------------------------------------------------
 Compute damage. Called from AssembleStress().
 ------------------------------------------------------------------------- */
void PairTlsph::ComputeDamage(const int i, const Matrix3d strain, const Matrix3d stress, Matrix3d &stress_damaged) {
	double *eff_plastic_strain = atom->eff_plastic_strain;
	double *eff_plastic_strain_rate = atom->eff_plastic_strain_rate;
	double *radius = atom->radius;
	double *damage = atom->damage;
	int *type = atom->type;
	int itype = type[i];
	double jc_failure_strain;
//double damage_gap, damage_rate;
	Matrix3d eye, stress_deviator;

	eye.setIdentity();
	stress_deviator = Deviator(stress);
	double pressure = -stress.trace() / 3.0;

	if (failureModel[itype].failure_max_principal_stress) {
		error->one(FLERR, "not yet implemented");
		/*
		 * maximum stress failure criterion:
		 */
		IsotropicMaxStressDamage(stress, Lookup[FAILURE_MAX_PRINCIPAL_STRESS_THRESHOLD][itype]);
	} else if (failureModel[itype].failure_max_principal_strain) {
		error->one(FLERR, "not yet implemented");
		/*
		 * maximum strain failure criterion:
		 */
		IsotropicMaxStrainDamage(strain, Lookup[FAILURE_MAX_PRINCIPAL_STRAIN_THRESHOLD][itype]);
	} else if (failureModel[itype].failure_max_plastic_strain) {
		if (eff_plastic_strain[i] >= Lookup[FAILURE_MAX_PLASTIC_STRAIN_THRESHOLD][itype]) {
			damage[i] = 1.0;
			//double damage_gap = 0.5 * Lookup[FAILURE_MAX_PLASTIC_STRAIN_THRESHOLD][itype];
			//damage[i] = (eff_plastic_strain[i] - Lookup[FAILURE_MAX_PLASTIC_STRAIN_THRESHOLD][itype]) / damage_gap;
		}
	} else if (failureModel[itype].failure_johnson_cook) {

		//cout << "this is stress deviator" << stress_deviator << endl;

		jc_failure_strain = JohnsonCookFailureStrain(pressure, stress_deviator, Lookup[FAILURE_JC_D1][itype],
				Lookup[FAILURE_JC_D2][itype], Lookup[FAILURE_JC_D3][itype], Lookup[FAILURE_JC_D4][itype],
				Lookup[FAILURE_JC_EPDOT0][itype], eff_plastic_strain_rate[i]);

		//cout << "plastic strain increment is " << plastic_strain_increment << "  jc fs is " << jc_failure_strain << endl;
		//printf("JC failure strain is: %f\n", jc_failure_strain);

		if (eff_plastic_strain[i] >= jc_failure_strain) {
			double damage_rate = Lookup[SIGNAL_VELOCITY][itype] / (100.0 * radius[i]);
			damage[i] += damage_rate * update->dt;
			//damage[i] = 1.0;
		}
	}

	/*
	 * Apply damage to integration point
	 */

//	damage[i] = MIN(damage[i], 0.8);
//
//	if (pressure > 0.0) { // compression: particle can carry compressive load but reduced shear
//		stress_damaged = -pressure * eye + (1.0 - damage[i]) * Deviator(stress);
//	} else { // tension: particle has reduced tensile and shear load bearing capability
//		stress_damaged = (1.0 - damage[i]) * (-pressure * eye + Deviator(stress));
//	}

}