/* ----------------------------------------------------------------------
 *
 *                    *** 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;
}