/* ----------------------------------------------------------------------
   LIGGGHTS - LAMMPS Improved for General Granular and Granular Heat
   Transfer Simulations

   LIGGGHTS is part of the CFDEMproject
   www.liggghts.com | www.cfdem.com

   Christoph Kloss, christoph.kloss@...1838...
   Copyright 2009-2012 JKU Linz
   Copyright 2012-     DCS Computing GmbH, Linz

   LIGGGHTS is based on LAMMPS
   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
   http://lammps.sandia.gov, Sandia National Laboratories
   Steve Plimpton, sjplimp@sandia.gov

   This software is distributed under the GNU General Public License.

   See the README file in the top-level directory.
------------------------------------------------------------------------- */

#include "string.h"
#include "stdlib.h"
#include "atom.h"
#include "update.h"
#include "respa.h"
#include "error.h"
#include "memory.h"
#include "modify.h"
#include "comm.h"
#include "math.h"
#include "vector_liggghts.h"
#include "mpi_liggghts.h"
#include "fix_Pforce.h"
#include "fix_property_atom.h"

using namespace LAMMPS_NS;
using namespace FixConst;

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

FixPforce::FixPforce(LAMMPS *lmp, int narg, char **arg) : Fix(lmp,narg,arg)
{
    fix_Pforce_ = NULL;


  peratom_flag = 1;
  size_peratom_cols = 3;
    vector_flag = 1;
    size_vector = 3;
    global_freq = 1;
    extvector = 1;
}

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

FixPforce::~FixPforce()
{

}

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

void FixPforce::post_create()
{
    // register dragforce
    if(!fix_Pforce_)
    {
        char* fixarg[11];
        fixarg[0]="Pforce";
        fixarg[1]="all";
        fixarg[2]="property/atom";
        fixarg[3]="Pforce";
        fixarg[4]="vector"; // 1 vector per particle to be registered
        fixarg[5]="yes";    // restart
        fixarg[6]="no";     // communicate ghost
        fixarg[7]="no";     // communicate rev
        fixarg[8]="0.001";
        fixarg[9]="0.";
        fixarg[10]="0.";
        fix_Pforce_ = modify->add_fix_property_atom(11,fixarg,style);
    }

}

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

void FixPforce::pre_delete(bool unfixflag)
{
    if(fix_Pforce_) modify->delete_fix("Pforce");
}

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

int FixPforce::setmask()
{
  int mask = 0;
  mask |= POST_FORCE;
  return mask;
}

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

void FixPforce::init()
{
    // make sure there is only one fix of this style
    if(modify->n_fixes_style(style) != 1)
      error->fix_error(FLERR,this,"More than one fix of style Pforce is not allowed");

    vectorZeroize3D(Pforce_total);
}

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

void FixPforce::post_force(int vflag)
{
  double **x = atom->x;
  double **f = atom->f;
  int *mask = atom->mask;
  int nlocal = atom->nlocal;
  double **Pforce = fix_Pforce_->array_atom;

  vectorZeroize3D(Pforce_total);

  // add Pforce to force vector
  for (int i = 0; i < nlocal; i++)
  {
    if (mask[i] & groupbit)
    {
        vectorAdd3D(f[i],Pforce[i],f[i]);
        vectorAdd3D(Pforce_total,Pforce[i],Pforce_total);
    }
  }
}

/* ----------------------------------------------------------------------
   return components of total force on fix group
------------------------------------------------------------------------- */

double FixPforce::compute_vector(int n)
{
 MPI_Sum_Vector(Pforce_total,3,world);
  return Pforce_total[n];
}