These are accelerated compute styles, which can be used if LAMMPS is diff --git a/doc/Section_commands.txt b/doc/Section_commands.txt index 284a097b..d341e093 100644 --- a/doc/Section_commands.txt +++ b/doc/Section_commands.txt @@ -485,6 +485,7 @@ each style or click on the style itself for a full description: "cna/atom"_compute_cna_atom.html, "com"_compute_com.html, "com/molecule"_compute_com_molecule.html, +"contact/atom"_compute_contact_atom.html, "coord/atom"_compute_coord_atom.html, "damage/atom"_compute_damage_atom.html, "dihedral/local"_compute_dihedral_local.html, diff --git a/doc/compute_contact_atom.html b/doc/compute_contact_atom.html new file mode 100644 index 00000000..6c35adf5 --- /dev/null +++ b/doc/compute_contact_atom.html @@ -0,0 +1,59 @@ + +
Syntax: +
+compute ID group-ID contact/atom ++
Examples: +
+compute 1 all contact/atom ++
Description: +
+Define a computation that calculates the number of contacts +for each atom in a group. +
+The contact number is defined for finite-size spherical particles as +the number of neighbor atoms which overlap the central particle, +meaning that their distance of separation is less than or equal to the +sum of the radii of the two particles. +
+The value of the contact number will be 0.0 for atoms not in the +specified compute group. +
+Output info: +
+This compute calculates a per-atom vector, whose values can be +accessed by any command that uses per-atom values from a compute as +input. See Section_howto 15 for an +overview of LAMMPS output options. +
+The per-atom vector values will be a number >= 0.0, as explained +above. +
+Restrictions: +
+This compute requires that atoms store a radius as defined by the +atom_style sphere command. +
+Related commands: +
+ +Default: none +
+ diff --git a/doc/compute_contact_atom.txt b/doc/compute_contact_atom.txt new file mode 100644 index 00000000..0b8b31a2 --- /dev/null +++ b/doc/compute_contact_atom.txt @@ -0,0 +1,54 @@ +"LAMMPS WWW Site"_lws - "LAMMPS Documentation"_ld - "LAMMPS Commands"_lc :c + +:link(lws,http://lammps.sandia.gov) +:link(ld,Manual.html) +:link(lc,Section_commands.html#comm) + +:line + +compute contact/atom command :h3 + +[Syntax:] + +compute ID group-ID contact/atom :pre + +ID, group-ID are documented in "compute"_compute.html command +contact/atom = style name of this compute command :ul + +[Examples:] + +compute 1 all contact/atom :pre + +[Description:] + +Define a computation that calculates the number of contacts +for each atom in a group. + +The contact number is defined for finite-size spherical particles as +the number of neighbor atoms which overlap the central particle, +meaning that their distance of separation is less than or equal to the +sum of the radii of the two particles. + +The value of the contact number will be 0.0 for atoms not in the +specified compute group. + +[Output info:] + +This compute calculates a per-atom vector, whose values can be +accessed by any command that uses per-atom values from a compute as +input. See "Section_howto 15"_Section_howto.html#howto_15 for an +overview of LAMMPS output options. + +The per-atom vector values will be a number >= 0.0, as explained +above. + +[Restrictions:] + +This compute requires that atoms store a radius as defined by the +"atom_style sphere"_atom_style.html command. + +[Related commands:] + +"compute coord/atom"_compute_coord_atom.html + +[Default:] none diff --git a/doc/dump.html b/doc/dump.html index 89e23a29..ca47b45c 100644 --- a/doc/dump.html +++ b/doc/dump.html @@ -38,7 +38,9 @@image args = discussed on dump image doc page
mesh/stl args = 'local' or 'ghost' or 'all'
- mesh/vtk args = one or more dump-identifiers
+ mesh/vtk args = zero or more keyword/ value pairs and one or more dump-identifiers
+ keywords = output
+ output values = face or interpolate
dump-identifier = 'stress' or 'id' or 'wear' or 'vel' or 'stresscomponents' or 'owner' or 'area'
decomposition/vtk args = none
@@ -348,6 +350,8 @@ all active meshes are dumped if you do not supply the optional list of
mesh IDs. By specifying list of mesh IDs you can explicitly choose
which meshes to dump. The group-ID is ignored. Again, a series of files
can be post-processed in Paraview , www.paraview.org
+Keyword output controlls if the data is written in a per-face manner
+or as interpolated values to VTK.
The decomposition/vtk style dumps the processor grid decomposition into a series of VTK files. No further args are expected. diff --git a/doc/dump.txt b/doc/dump.txt index 9ef9ad55..2eecdcdc 100644 --- a/doc/dump.txt +++ b/doc/dump.txt @@ -28,7 +28,9 @@ args = list of arguments for a particular style :l {image} args = discussed on "dump image"_dump_image.html doc page :pre {mesh/stl} args = 'local' or 'ghost' or 'all' - {mesh/vtk} args = one or more dump-identifiers + {mesh/vtk} args = zero or more keyword/ value pairs and one or more dump-identifiers + keywords = {output} + {output} values = face or interpolate dump-identifier = 'stress' or 'id' or 'wear' or 'vel' or 'stresscomponents' or 'owner' or 'area' {decomposition/vtk} args = none :pre @@ -337,6 +339,8 @@ all active meshes are dumped if you do not supply the optional list of mesh IDs. By specifying list of mesh IDs you can explicitly choose which meshes to dump. The group-ID is ignored. Again, a series of files can be post-processed in Paraview , www.paraview.org +Keyword {output} controlls if the data is written in a per-face manner +or as interpolated values to VTK. The {decomposition/vtk} style dumps the processor grid decomposition into a series of VTK files. No further args are expected. diff --git a/doc/fix_heat_gran.html b/doc/fix_heat_gran.html index bd4f62ca..d7b41b35 100644 --- a/doc/fix_heat_gran.html +++ b/doc/fix_heat_gran.html @@ -29,7 +29,7 @@
area_correction values = on or off +area_correction values = yes or noExamples:
diff --git a/doc/fix_heat_gran.txt b/doc/fix_heat_gran.txt index f6947474..a984bf36 100644 --- a/doc/fix_heat_gran.txt +++ b/doc/fix_heat_gran.txt @@ -19,7 +19,7 @@ initial_temperature = obligatory keyword :l T0 = initial (default) temperature for the particles :l zero or more keyword/value pairs may be appended :l keyword = {area_correction} :l - {area_correction} values = {on} or {off} :pre + {area_correction} values = {yes} or {no} :pre [Examples:] diff --git a/doc/fix_insert_pack.html b/doc/fix_insert_pack.html index bbbc279a..a376f2d7 100644 --- a/doc/fix_insert_pack.html +++ b/doc/fix_insert_pack.html @@ -31,9 +31,10 @@
maxattempt value = ma
+ verbose = yes or no
+ maxattempt value = ma
ma = max # of insertion attempts per atom (positive integer)
insert_every value = once or ie
once = value to signalise that isertion takes place only once (the step after the fix is issued)
@@ -79,6 +80,9 @@
Insert particles into a granular run either once or every few timesteps
within the specified region, as defined via the region keyword.
+The verbose keyword controlls whether statistics about particle
+insertion is output to the screen each time particles are inserted.
+
This command must use the distributiontemplate keyword to refer to a
fix_particledistribution_discrete (defined by dist-fix-ID) that defines the properties of the
inserted particles.
diff --git a/doc/fix_insert_pack.txt b/doc/fix_insert_pack.txt
index 6ca32326..c536439c 100644
--- a/doc/fix_insert_pack.txt
+++ b/doc/fix_insert_pack.txt
@@ -20,7 +20,8 @@ seed_value = random # seed (positive integer) :l
distributiontemplate = obligatory keyword :l
dist-ID = ID of a "fix_particledistribution_discrete"_fix_particledistribution_discrete.html to be used for particle insertion :l
one or more general keyword/value pairs can be appended :l
-general_keywords = {maxattampt} or {insert_every} or {overlapcheck} or {all_in} or {random_distribute} or {vel} or {omega} :l
+general_keywords = {verbose} or {maxattampt} or {insert_every} or {overlapcheck} or {all_in} or {random_distribute} or {vel} or {omega} :l
+ {verbose} = yes or no
{maxattempt} value = ma
ma = max # of insertion attempts per atom (positive integer)
{insert_every} value = once or ie
@@ -64,6 +65,9 @@ fix ins all insert/pack seed 1001 distributiontemplate pdd1 insert_every once ov
Insert particles into a granular run either once or every few timesteps
within the specified region, as defined via the {region} keyword.
+The {verbose} keyword controlls whether statistics about particle
+insertion is output to the screen each time particles are inserted.
+
This command must use the distributiontemplate keyword to refer to a
"fix_particledistribution_discrete"_fix_particledistribution_discrete.html (defined by dist-fix-ID) that defines the properties of the
inserted particles.
diff --git a/doc/fix_insert_rate_region.html b/doc/fix_insert_rate_region.html
index 944b280e..1134681c 100644
--- a/doc/fix_insert_rate_region.html
+++ b/doc/fix_insert_rate_region.html
@@ -31,9 +31,10 @@
maxattempt value = ma
+ verbose = yes or no
+ maxattempt value = ma
ma = max # of insertion attempts per atom (positive integer)
nparticles values = np
np = number of particles to insert (positive integer)
@@ -80,6 +81,9 @@
Insert particles into a granular run every few timesteps
within a specified region at a specified rate.
+The verbose keyword controlls whether statistics about particle
+insertion is output to the screen each time particles are inserted.
+
This command must use the region keyword to define an insertion
volume. The specified region must have been previously defined with a
region command. Dynamic regions are not supported as
diff --git a/doc/fix_insert_rate_region.txt b/doc/fix_insert_rate_region.txt
index 4e53e000..c452e793 100644
--- a/doc/fix_insert_rate_region.txt
+++ b/doc/fix_insert_rate_region.txt
@@ -20,7 +20,8 @@ seed_value = random # seed (positive integer) :l
distributiontemplate = obligatory keyword :l
dist-ID = ID of a "fix_particledistribution_discrete"_fix_particledistribution_discrete.html to be used for particle insertion :l
one or more general keyword/value pairs can be appended :l
-general_keywords = {maxattampt} or {nparticles} or {mass} or {particlerate} or {massrate} or {insert_every} or {overlapcheck} or {all_in} or {random_distribute} or {vel} or {omega} :l
+general_keywords = {verbose} or {maxattampt} or {nparticles} or {mass} or {particlerate} or {massrate} or {insert_every} or {overlapcheck} or {all_in} or {random_distribute} or {vel} or {omega} :l
+ {verbose} = yes or no
{maxattempt} value = ma
ma = max # of insertion attempts per atom (positive integer)
{nparticles} values = np
@@ -66,6 +67,9 @@ fix ins all insert/rate/region seed 1001 distributiontemplate pdd1 nparticles 10
Insert particles into a granular run every few timesteps
within a specified region at a specified rate.
+The {verbose} keyword controlls whether statistics about particle
+insertion is output to the screen each time particles are inserted.
+
This command must use the {region} keyword to define an insertion
volume. The specified region must have been previously defined with a
"region"_region.html command. Dynamic regions are not supported as
diff --git a/doc/fix_insert_stream.html b/doc/fix_insert_stream.html
index f586ac57..9bf3bc2a 100644
--- a/doc/fix_insert_stream.html
+++ b/doc/fix_insert_stream.html
@@ -31,9 +31,10 @@
maxattempt value = ma
+ verbose = yes or no
+ maxattempt value = ma
ma = max # of insertion attempts per atom (positive integer)
nparticles values = np
np = number of particles to insert (positive integer)
@@ -82,6 +83,9 @@
within a specified region until either np particles have been inserted
or the desired particle mass (mp) has been reached.
+The verbose keyword controlls whether statistics about particle
+insertion is output to the screen each time particles are inserted.
+
Each timestep particles are inserted, they are placed randomly inside
the insertion volume so as to mimic a stream of poured particles. The
insertion volume is generated by extruding the insertion face as
diff --git a/doc/fix_insert_stream.txt b/doc/fix_insert_stream.txt
index 1a8cd8b8..2a2c9a43 100644
--- a/doc/fix_insert_stream.txt
+++ b/doc/fix_insert_stream.txt
@@ -20,7 +20,8 @@ seed_value = random # seed (positive integer) :l
distributiontemplate = obligatory keyword :l
dist-ID = ID of a "fix_particledistribution_discrete"_fix_particledistribution_discrete.html to be used for particle insertion :l
one or more general keyword/value pairs can be appended :l
-general_keywords = {maxattampt} or {nparticles} or {mass} or {particlerate} or {massrate} or {insert_every} or {overlapcheck} or {all_in} or {random_distribute} or {vel} or {omega} :l
+general_keywords = {verbose} or {maxattampt} or {nparticles} or {mass} or {particlerate} or {massrate} or {insert_every} or {overlapcheck} or {all_in} or {random_distribute} or {vel} or {omega} :l
+ {verbose} = yes or no
{maxattempt} value = ma
ma = max # of insertion attempts per atom (positive integer)
{nparticles} values = np
@@ -67,6 +68,9 @@ Insert particles into a granular run either once or every few timesteps
within a specified region until either np particles have been inserted
or the desired particle mass (mp) has been reached.
+The {verbose} keyword controlls whether statistics about particle
+insertion is output to the screen each time particles are inserted.
+
Each timestep particles are inserted, they are placed randomly inside
the insertion volume so as to mimic a stream of poured particles. The
insertion volume is generated by extruding the insertion face as
diff --git a/src/compute_contact_atom.cpp b/src/compute_contact_atom.cpp
new file mode 100644
index 00000000..ef3d4ff5
--- /dev/null
+++ b/src/compute_contact_atom.cpp
@@ -0,0 +1,197 @@
+/* ----------------------------------------------------------------------
+ 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 "math.h"
+#include "string.h"
+#include "stdlib.h"
+#include "compute_contact_atom.h"
+#include "atom.h"
+#include "update.h"
+#include "modify.h"
+#include "neighbor.h"
+#include "neigh_list.h"
+#include "neigh_request.h"
+#include "force.h"
+#include "pair.h"
+#include "comm.h"
+#include "memory.h"
+#include "error.h"
+
+using namespace LAMMPS_NS;
+
+/* ---------------------------------------------------------------------- */
+
+ComputeContactAtom::ComputeContactAtom(LAMMPS *lmp, int narg, char **arg) :
+ Compute(lmp, narg, arg)
+{
+ if (narg != 3) error->all(FLERR,"Illegal compute contact/atom command");
+
+ peratom_flag = 1;
+ size_peratom_cols = 0;
+ comm_reverse = 1;
+
+ nmax = 0;
+ contact = NULL;
+
+ // error checks
+
+ if (!atom->sphere_flag)
+ error->all(FLERR,"Compute contact/atom requires atom style sphere");
+}
+
+/* ---------------------------------------------------------------------- */
+
+ComputeContactAtom::~ComputeContactAtom()
+{
+ memory->destroy(contact);
+}
+
+/* ---------------------------------------------------------------------- */
+
+void ComputeContactAtom::init()
+{
+ if (force->pair == NULL)
+ error->all(FLERR,"Compute contact/atom requires a pair style be defined");
+
+ int count = 0;
+ for (int i = 0; i < modify->ncompute; i++)
+ if (strcmp(modify->compute[i]->style,"contact/atom") == 0) count++;
+ if (count > 1 && comm->me == 0)
+ error->warning(FLERR,"More than one compute contact/atom");
+
+ // need an occasional neighbor list
+
+ int irequest = neighbor->request((void *) this);
+ neighbor->requests[irequest]->half = 0;
+ neighbor->requests[irequest]->gran = 1;
+ neighbor->requests[irequest]->pair = 0;
+ neighbor->requests[irequest]->compute = 1;
+ neighbor->requests[irequest]->occasional = 1;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void ComputeContactAtom::init_list(int id, NeighList *ptr)
+{
+ list = ptr;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void ComputeContactAtom::compute_peratom()
+{
+ int i,j,ii,jj,inum,jnum;
+ double xtmp,ytmp,ztmp,delx,dely,delz,rsq;
+ double radi,radsum,radsumsq;
+ int *ilist,*jlist,*numneigh,**firstneigh;
+
+ invoked_peratom = update->ntimestep;
+
+ // grow contact array if necessary
+
+ if (atom->nmax > nmax) {
+ memory->destroy(contact);
+ nmax = atom->nmax;
+ memory->create(contact,nmax,"contact/atom:contact");
+ vector_atom = contact;
+ }
+
+ // invoke neighbor list (will copy or build if necessary)
+
+ neighbor->build_one(list->index);
+
+ inum = list->inum;
+ ilist = list->ilist;
+ numneigh = list->numneigh;
+ firstneigh = list->firstneigh;
+
+ // compute number of contacts for each atom in group
+ // contact if distance <= sum of radii
+ // tally for both I and J
+
+ double **x = atom->x;
+ double *radius = atom->radius;
+ int *mask = atom->mask;
+ int nlocal = atom->nlocal;
+ int nall = nlocal + atom->nghost;
+
+ for (i = 0; i < nall; i++) contact[i] = 0.0;
+
+ for (ii = 0; ii < inum; ii++) {
+ i = ilist[ii];
+ if (mask[i] & groupbit) {
+ xtmp = x[i][0];
+ ytmp = x[i][1];
+ ztmp = x[i][2];
+ radi = radius[i];
+ jlist = firstneigh[i];
+ jnum = numneigh[i];
+
+ for (jj = 0; jj < jnum; jj++) {
+ j = jlist[jj];
+ j &= NEIGHMASK;
+
+ delx = xtmp - x[j][0];
+ dely = ytmp - x[j][1];
+ delz = ztmp - x[j][2];
+ rsq = delx*delx + dely*dely + delz*delz;
+ radsum = radi + radius[j];
+ radsumsq = radsum*radsum;
+ if (rsq <= radsumsq) {
+ contact[i] += 1.0;
+ contact[j] += 1.0;
+ }
+ }
+ }
+ }
+
+ // communicate ghost atom counts between neighbor procs if necessary
+
+ if (force->newton_pair) comm->reverse_comm_compute(this);
+}
+
+/* ---------------------------------------------------------------------- */
+
+int ComputeContactAtom::pack_reverse_comm(int n, int first, double *buf)
+{
+ int i,m,last;
+
+ m = 0;
+ last = first + n;
+ for (i = first; i < last; i++)
+ buf[m++] = contact[i];
+ return 1;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void ComputeContactAtom::unpack_reverse_comm(int n, int *list, double *buf)
+{
+ int i,j,m;
+
+ m = 0;
+ for (i = 0; i < n; i++) {
+ j = list[i];
+ contact[j] += buf[m++];
+ }
+}
+
+/* ----------------------------------------------------------------------
+ memory usage of local atom-based array
+------------------------------------------------------------------------- */
+
+double ComputeContactAtom::memory_usage()
+{
+ double bytes = nmax * sizeof(double);
+ return bytes;
+}
diff --git a/src/compute_contact_atom.h b/src/compute_contact_atom.h
new file mode 100644
index 00000000..3747d42c
--- /dev/null
+++ b/src/compute_contact_atom.h
@@ -0,0 +1,66 @@
+/* ----------------------------------------------------------------------
+ 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.
+------------------------------------------------------------------------- */
+
+#ifdef COMPUTE_CLASS
+
+ComputeStyle(coord/gran,ComputeContactAtom)
+ComputeStyle(contact/atom,ComputeContactAtom)
+
+#else
+
+#ifndef LMP_COMPUTE_CONTACT_ATOM_H
+#define LMP_COMPUTE_CONTACT_ATOM_H
+
+#include "compute.h"
+
+namespace LAMMPS_NS {
+
+class ComputeContactAtom : public Compute {
+ public:
+ ComputeContactAtom(class LAMMPS *, int, char **);
+ ~ComputeContactAtom();
+ void init();
+ void init_list(int, class NeighList *);
+ void compute_peratom();
+ int pack_reverse_comm(int, int, double *);
+ void unpack_reverse_comm(int, int *, double *);
+ double memory_usage();
+
+ private:
+ int nmax;
+ class NeighList *list;
+ double *contact;
+};
+
+}
+
+#endif
+#endif
+
+/* ERROR/WARNING messages:
+
+E: Illegal ... command
+
+Self-explanatory. Check the input script syntax and compare to the
+documentation for the command. You can use -echo screen as a
+command-line option when running LAMMPS to see the offending line.
+
+E: Compute contact/atom requires a pair style be defined
+
+Self-explantory.
+
+W: More than one compute contact/atom
+
+It is not efficient to use compute contact/atom more than once.
+
+*/
diff --git a/src/compute_pair_gran_local.cpp b/src/compute_pair_gran_local.cpp
index bd7bbd97..7070bc73 100644
--- a/src/compute_pair_gran_local.cpp
+++ b/src/compute_pair_gran_local.cpp
@@ -452,7 +452,7 @@ void ComputePairGranLocal::add_heat(int i,int j,double hf)
add data from particle-wall contact on this proc
------------------------------------------------------------------------- */
-void ComputePairGranLocal::add_wall_1(int iFMG,int iTri,int iP,double *contact_point)
+void ComputePairGranLocal::add_wall_1(int iFMG,int idTri,int iP,double *contact_point)
{
if (!(atom->mask[iP] & groupbit)) return;
@@ -470,7 +470,7 @@ void ComputePairGranLocal::add_wall_1(int iFMG,int iTri,int iP,double *contact_p
if(idflag)
{
array[ipair][n++] = static_cast(iFMG);
- array[ipair][n++] = static_cast(iTri);
+ array[ipair][n++] = static_cast(idTri);
}
}
diff --git a/src/dump_custom.cpp b/src/dump_custom.cpp
index aab2d960..fd596c50 100644
--- a/src/dump_custom.cpp
+++ b/src/dump_custom.cpp
@@ -1365,6 +1365,15 @@ int DumpCustom::modify_param(int narg, char **arg)
return 2;
}
+ if (strcmp(arg[0],"label") == 0) {
+ if (narg < 2) error->all(FLERR,"Illegal dump_modify command [label]");
+ delete [] label;
+ int n = strlen(arg[1]) + 1;
+ label = new char[n];
+ strcpy(label,arg[1]);
+ return 2;
+ }
+
if (strcmp(arg[0],"element") == 0) {
if (narg < ntypes+1)
error->all(FLERR,"Dump modify element names do not match atom types");
diff --git a/src/dump_euler_vtk.cpp b/src/dump_euler_vtk.cpp
new file mode 100644
index 00000000..448c4294
--- /dev/null
+++ b/src/dump_euler_vtk.cpp
@@ -0,0 +1,208 @@
+/* ----------------------------------------------------------------------
+ LIGGGHTS - LAMMPS Improved for General Granular and Granular Heat
+ Transfer Simulations
+
+ LIGGGHTS is part of the CFDEMproject
+ www.liggghts.com | www.cfdem.com
+
+ Christoph Kloss, christoph.kloss@cfdem.com
+ Copyright 2009-2012 JKU Linz
+ Copyright 2012- DCS Computing GmbH, Linz
+
+ LIGGGHTS is based on LAMMPS
+ LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+ http://lammps.sandia.gov, Sandia National Laboratories
+ Steve Plimpton, sjplimp@sandia.gov
+
+ This software is distributed under the GNU General Public License.
+
+ See the README file in the top-level directory.
+------------------------------------------------------------------------- */
+
+#include "string.h"
+#include "dump_euler_vtk.h"
+#include "fix_ave_euler.h"
+#include "domain.h"
+#include "atom.h"
+#include "update.h"
+#include "group.h"
+#include "error.h"
+#include "fix.h"
+#include "modify.h"
+#include "comm.h"
+#include
+
+using namespace LAMMPS_NS;
+
+/* ---------------------------------------------------------------------- */
+
+DumpEulerVTK::DumpEulerVTK(LAMMPS *lmp, int narg, char **arg) : Dump(lmp, narg, arg),
+ fix_euler_(0)
+{
+ if (narg < 5)
+ error->all(FLERR,"Illegal dump pic/vtk command");
+
+ // CURRENTLY ONLY PROC 0 writes
+
+ format_default = NULL;
+}
+
+/* ---------------------------------------------------------------------- */
+
+DumpEulerVTK::~DumpEulerVTK()
+{
+}
+
+/* ---------------------------------------------------------------------- */
+
+void DumpEulerVTK::init_style()
+{
+ fix_euler_ = static_cast(modify->find_fix_style("ave/euler",0));
+ if(!fix_euler_)
+ error->all(FLERR,"Illegal dump pic/vtk command, need a fix ave/euler");
+
+ // multifile=1; // 0 = one big file, 1 = one file per timestep
+ // multiproc=0; // 0 = proc 0 writes for all, 1 = one file/proc
+ if (multifile != 1)
+ error->all(FLERR,"You should use a filename like 'dump*.vtk' for the 'dump pic/vtk' command to produce one file per time-step");
+ if (multiproc != 0)
+ error->all(FLERR,"Your 'dump pic/vtk' command is writing one file per processor, where all the files contain the same data");
+
+// if (domain->triclinic == 1)
+// error->all(FLERR,"Can not dump VTK files for triclinic box");
+ if (binary)
+ error->all(FLERR,"Can not dump VTK files in binary mode");
+
+ // node center (3), av vel (3), volume fraction, stress
+ size_one = 8;
+
+ delete [] format;
+}
+
+/* ---------------------------------------------------------------------- */
+
+int DumpEulerVTK::modify_param(int narg, char **arg)
+{
+ error->warning(FLERR,"dump_modify keyword is not supported by 'dump pic/vtk' and is thus ignored");
+ return 0;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void DumpEulerVTK::write_header(bigint ndump)
+{
+ write_header_ascii(ndump);
+}
+
+void DumpEulerVTK::write_header_ascii(bigint ndump)
+{
+ if (comm->me!=0) return;
+ fprintf(fp,"# vtk DataFile Version 2.0\nLIGGGHTS mesh/VTK export\nASCII\n");
+}
+
+/* ---------------------------------------------------------------------- */
+//NP number of tris in the dump file
+int DumpEulerVTK::count()
+{
+ return fix_euler_->ncells();
+}
+
+/* ---------------------------------------------------------------------- */
+
+void DumpEulerVTK::pack(int *ids)
+{
+ int m = 0;
+
+ // have to stick with this order (all per-element props)
+ // as multiple procs pack
+
+ int ncells = fix_euler_->ncells();
+
+ for(int i = 0; i < ncells; i++)
+ {
+ buf[m++] = fix_euler_->cell_center(i,0);
+ buf[m++] = fix_euler_->cell_center(i,1);
+ buf[m++] = fix_euler_->cell_center(i,2);
+
+ buf[m++] = fix_euler_->cell_v_av(i,0);
+ buf[m++] = fix_euler_->cell_v_av(i,1);
+ buf[m++] = fix_euler_->cell_v_av(i,2);
+
+ buf[m++] = fix_euler_->cell_vol_fr(i);
+
+ buf[m++] = fix_euler_->cell_pressure(i);
+ }
+ return ;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void DumpEulerVTK::write_data(int n, double *mybuf)
+{
+ write_data_ascii(n,mybuf);
+}
+
+void DumpEulerVTK::write_data_ascii(int n, double *mybuf)
+{
+ //only proc 0 writes
+ if (comm->me != 0) return;
+
+ int m, buf_pos;
+
+ // n is the number of elements
+
+ /*NL*///fprintf(screen,"WRITING ITEM at step %d proc %d with n %d\n",update->ntimestep,comm->me,n);
+
+ // write point data
+ fprintf(fp,"DATASET POLYDATA\nPOINTS %d float\n",n);
+ m = 0;
+ buf_pos = 0;
+ for (int i = 0; i < n; i++)
+ {
+ fprintf(fp,"%f %f %f\n",mybuf[m],mybuf[m+1],mybuf[m+2]);
+ m += size_one ;
+ }
+ buf_pos += 3;
+
+ // write polygon data
+ fprintf(fp,"VERTICES %d %d\n",n,2*n);
+ for (int i = 0; i < n; i++)
+ {
+ fprintf(fp,"%d %d\n",1,i);
+ }
+
+ // write point data header
+ fprintf(fp,"POINT_DATA %d\n",n);
+
+ // write cell data
+
+ fprintf(fp,"VECTORS v_avg float\n");
+ m = buf_pos;
+ for (int i = 0; i < n; i++)
+ {
+ fprintf(fp,"%f %f %f\n",mybuf[m],mybuf[m+1],mybuf[m+2]);
+ m += size_one;
+ }
+ buf_pos += 3;
+
+ fprintf(fp,"SCALARS volumefraction float 1\nLOOKUP_TABLE default\n");
+ m = buf_pos;
+ for (int i = 0; i < n; i++)
+ {
+ fprintf(fp,"%f\n",mybuf[m]);
+ m += size_one;
+ }
+ buf_pos++;
+
+ fprintf(fp,"SCALARS pressure float 1\nLOOKUP_TABLE default\n");
+ m = buf_pos;
+ for (int i = 0; i < n; i++)
+ {
+ fprintf(fp,"%f\n",mybuf[m]);
+ m += size_one;
+ }
+ buf_pos++;
+
+ // footer not needed
+ // if would be needed, would do like in dump stl
+}
diff --git a/src/dump_euler_vtk.h b/src/dump_euler_vtk.h
new file mode 100644
index 00000000..e4c8fa7f
--- /dev/null
+++ b/src/dump_euler_vtk.h
@@ -0,0 +1,60 @@
+/* ----------------------------------------------------------------------
+ LIGGGHTS - LAMMPS Improved for General Granular and Granular Heat
+ Transfer Simulations
+
+ LIGGGHTS is part of the CFDEMproject
+ www.liggghts.com | www.cfdem.com
+
+ Christoph Kloss, christoph.kloss@cfdem.com
+ Copyright 2009-2012 JKU Linz
+ Copyright 2012- DCS Computing GmbH, Linz
+
+ LIGGGHTS is based on LAMMPS
+ LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+ http://lammps.sandia.gov, Sandia National Laboratories
+ Steve Plimpton, sjplimp@sandia.gov
+
+ This software is distributed under the GNU General Public License.
+
+ See the README file in the top-level directory.
+------------------------------------------------------------------------- */
+
+#ifdef DUMP_CLASS
+
+DumpStyle(euler/vtk,DumpEulerVTK)
+
+#else
+
+#ifndef LMP_DUMP_EULER_VTK_H
+#define LMP_DUMP_EULER_VTK_H
+
+#include "dump.h"
+
+namespace LAMMPS_NS {
+
+class DumpEulerVTK : public Dump {
+
+ public:
+
+ DumpEulerVTK(LAMMPS *, int, char**);
+ virtual ~DumpEulerVTK();
+ void init_style();
+
+ private:
+
+ class FixAveEuler *fix_euler_;
+
+ int modify_param(int, char **);
+ void write_header(bigint ndump);
+ int count();
+ void pack(int *);
+ void write_data(int, double *);
+
+ void write_header_ascii(bigint ndump);
+ void write_data_ascii(int n, double *mybuf);
+};
+
+}
+
+#endif
+#endif
diff --git a/src/dump_mesh_vtk.cpp b/src/dump_mesh_vtk.cpp
index a0f2f4ed..378c586f 100644
--- a/src/dump_mesh_vtk.cpp
+++ b/src/dump_mesh_vtk.cpp
@@ -19,6 +19,11 @@
See the README file in the top-level directory.
------------------------------------------------------------------------- */
+/* ----------------------------------------------------------------------
+ Contributing author for interpolated output:
+ Felix Kleinfeldt (OVGU Magdeburg)
+------------------------------------------------------------------------- */
+
#include "string.h"
#include "dump_mesh_vtk.h"
#include "tri_mesh.h"
@@ -35,7 +40,8 @@
using namespace LAMMPS_NS;
-enum{
+enum
+{
DUMP_STRESS = 1,
DUMP_STRESSCOMPONENTS = 2,
DUMP_ID = 4,
@@ -44,8 +50,13 @@ enum{
DUMP_TEMP = 32,
DUMP_OWNER = 64,
DUMP_AREA = 128
- };
+};
+enum
+{
+ DUMP_POINTS = 1,
+ DUMP_FACE = 2
+};
/* ---------------------------------------------------------------------- */
@@ -93,12 +104,22 @@ DumpMeshVTK::DumpMeshVTK(LAMMPS *lmp, int narg, char **arg) : Dump(lmp, narg, ar
int iarg = 5;
dump_what_ = 0;
+ dataMode_;
bool hasargs = true;
while (iarg < narg && hasargs)
{
hasargs = false;
- if(strcmp(arg[iarg],"stress")==0)
+ if(strcmp(arg[iarg],"output")==0)
+ {
+ if (iarg+2 > narg) error->all(FLERR,"Dump mesh/vtk: not enough arguments for 'interpolate'");
+ if(strcmp(arg[iarg+1],"face")==0) dataMode_ = 0;
+ else if(strcmp(arg[iarg+1],"interpolate")==0) dataMode_ = 1;
+ else error->all(FLERR,"Dump mesh/vtk: wrong arrgument for 'interpolate'");
+ iarg += 2;
+ hasargs = true;
+ }
+ else if(strcmp(arg[iarg],"stress")==0)
{
dump_what_ |= DUMP_STRESS;
iarg++;
@@ -411,7 +432,240 @@ void DumpMeshVTK::write_data(int n, double *mybuf)
write_data_ascii(n_all_/size_one,buf_all_);
}
+/* ---------------------------------------------------------------------- */
+
void DumpMeshVTK::write_data_ascii(int n, double *mybuf)
+{
+ if(dataMode_ == 0)
+ write_data_ascii_face(n,mybuf);
+ else if(dataMode_ == 1)
+ write_data_ascii_point(n,mybuf);
+ return;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void DumpMeshVTK::write_data_ascii_point(int n, double *mybuf)
+{
+ int m, buf_pos;
+
+ // n is the number of tri
+
+ m = 0;
+ buf_pos = 0;
+
+ ScalarContainer points;
+ ScalarContainer tri_points;
+ bool add;
+ for (int i = 0; i < n; i++)
+ {
+ for (int j=0;j<3;j++)
+ {
+ add = true;
+ for (int k=0; k < points.size(); k+=3)
+ {
+ if (points.get(k)==mybuf[m+0] && points.get(k+1)==mybuf[m+1] && points.get(k+2)==mybuf[m+2])
+ {
+ add=false;
+ tri_points.add(k/3);
+ break;
+ }
+ }
+ if (add)
+ {
+ for (int k=0;k<3;k++)
+ points.add(mybuf[m+k]);
+ int max=tri_points.max();
+ if (max<0) max=-1;
+ tri_points.add(max+1);
+ }
+ m+=3;
+ }
+ m += size_one-9;
+ }
+
+ buf_pos += 9;
+
+ // points_neightri
+ class ScalarContainer **points_neightri;
+ points_neightri = new ScalarContainer*[(int)points.size()/3];
+ for (int i=0; i < points.size()/3; i++)
+ points_neightri[i] = new ScalarContainer;
+ for (int i=0; i < 3*n; i+=3)
+ {
+ for (int j=0; j<3;j++)
+ points_neightri[tri_points.get(i+j)]->add(i/3);
+ }
+
+
+ // write point data
+ fprintf(fp,"DATASET UNSTRUCTURED_GRID\nPOINTS %d float\n", points.size()/3);
+ for (int i=0; i < points.size(); i+=3)
+ fprintf(fp,"%f %f %f\n",points.get(i+0),points.get(i+1),points.get(i+2));
+
+ // write polygon data
+ fprintf(fp,"CELLS %d %d\n",n,4*n);
+ for (int i = 0; i < 3*n; i+=3) fprintf(fp,"%d %d %d %d\n",3,tri_points.get(i+0),tri_points.get(i+1),tri_points.get(i+2));
+ // write cell types
+ fprintf(fp,"CELL_TYPES %d\n",n);
+ for (int i = 0; i < n; i++)
+ fprintf(fp,"5\n");
+
+ // write point data header
+ fprintf(fp,"POINT_DATA %d\n",points.size()/3);
+
+ // write cell data
+
+ if(dump_what_ & DUMP_STRESS)
+ {
+ // write pressure and shear stress
+ fprintf(fp,"SCALARS pressure float 1\nLOOKUP_TABLE default\n");
+ m = buf_pos;
+ for (int i = 0; i < points.size()/3; i++)
+ {
+ double helper=0;
+ for (int j=0; j < points_neightri[i]->size();j++) helper += mybuf[m + points_neightri[i]->get(j)*size_one];
+ helper /= points_neightri[i]->size();
+ fprintf(fp,"%f\n",helper);
+ }
+ buf_pos++;
+
+ // write shear stress
+ fprintf(fp,"SCALARS shearstress float 1\nLOOKUP_TABLE default\n");
+ m = buf_pos;
+ for (int i = 0; i < points.size()/3; i++)
+ {
+ double helper=0;
+ for (int j=0; j < points_neightri[i]->size();j++) helper += mybuf[m + points_neightri[i]->get(j)*size_one];
+ helper /= points_neightri[i]->size();
+ fprintf(fp,"%f\n",helper);
+ }
+ buf_pos++;
+ }
+ if(dump_what_ & DUMP_STRESSCOMPONENTS)
+ {
+ //write x y z stress component
+ fprintf(fp,"VECTORS stress float\n");
+ m = buf_pos;
+ for (int i = 0; i < points.size()/3; i++)
+ {
+ double helper1=0, helper2=0, helper3=0;
+ for (int j=0; j < points_neightri[i]->size();j++)
+ {
+ helper1 += mybuf[m + points_neightri[i]->get(j)*size_one];
+ helper2 += mybuf[m+1 + points_neightri[i]->get(j)*size_one];
+ helper3 += mybuf[m+2 + points_neightri[i]->get(j)*size_one];
+ }
+ helper1 /= points_neightri[i]->size();
+ helper2 /= points_neightri[i]->size();
+ helper3 /= points_neightri[i]->size();
+ fprintf(fp,"%f %f %f\n",helper1,helper2,helper3);
+ }
+ buf_pos += 3;
+ }
+
+ if(dump_what_ & DUMP_ID)
+ {
+ // write id
+ fprintf(fp,"SCALARS meshid float 1\nLOOKUP_TABLE default\n");
+ m = buf_pos;
+ for (int i = 0; i < points.size()/3; i++)
+ {
+ double helper=0;
+ for (int j=0; j < points_neightri[i]->size();j++) helper += mybuf[m + points_neightri[i]->get(j)*size_one];
+ helper /= points_neightri[i]->size();
+ fprintf(fp,"%f\n",helper);
+ }
+ buf_pos++;
+ }
+
+ if(dump_what_ & DUMP_VEL)
+ {
+ //write vel data
+ fprintf(fp,"VECTORS v float\n");
+ m = buf_pos;
+ for (int i = 0; i < points.size()/3; i++)
+ {
+ double helper1=0, helper2=0, helper3=0;
+ for (int j=0; j < points_neightri[i]->size();j++)
+ {
+ helper1 += mybuf[m + points_neightri[i]->get(j)*size_one];
+ helper2 += mybuf[m+1 + points_neightri[i]->get(j)*size_one];
+ helper3 += mybuf[m+2 + points_neightri[i]->get(j)*size_one];
+ }
+ helper1 /= points_neightri[i]->size();
+ helper2 /= points_neightri[i]->size();
+ helper3 /= points_neightri[i]->size();
+ fprintf(fp,"%f %f %f\n",helper1,helper2,helper3);
+ }
+ buf_pos += 3;
+ }
+
+ if(dump_what_ & DUMP_WEAR)
+ {
+ //write wear data
+ fprintf(fp,"SCALARS wear float 1\nLOOKUP_TABLE default\n");
+ m = buf_pos;
+ for (int i = 0; i < points.size()/3; i++)
+ {
+ double helper=0;
+ for (int j=0; j < points_neightri[i]->size();j++) helper += mybuf[m + points_neightri[i]->get(j)*size_one];
+ helper /= points_neightri[i]->size();
+ fprintf(fp,"%f\n",helper);
+ }
+ buf_pos++;
+ }
+
+ if(dump_what_ & DUMP_TEMP)
+ {
+ //write wear data
+ fprintf(fp,"SCALARS Temp float 1\nLOOKUP_TABLE default\n");
+ m = buf_pos;
+ for (int i = 0; i < points.size()/3; i++)
+ {
+ double helper=0;
+ for (int j=0; j < points_neightri[i]->size();j++) helper += mybuf[m + points_neightri[i]->get(j)*size_one];
+ helper /= points_neightri[i]->size();
+ fprintf(fp,"%f\n",helper);
+ }
+ buf_pos++;
+ }
+
+ if(dump_what_ & DUMP_OWNER)
+ {
+ //write owner data
+ fprintf(fp,"SCALARS owner float 1\nLOOKUP_TABLE default\n");
+ m = buf_pos;
+ for (int i = 0; i < points.size()/3; i++)
+ {
+ double helper=0;
+ for (int j=0; j < points_neightri[i]->size();j++) helper += mybuf[m + points_neightri[i]->get(j)*size_one];
+ helper /= points_neightri[i]->size();
+ fprintf(fp,"%f\n",helper);
+ }
+ buf_pos++;
+ }
+
+ if(dump_what_ & DUMP_AREA)
+ {
+ //write area data
+ fprintf(fp,"SCALARS area float 1\nLOOKUP_TABLE default\n");
+ m = buf_pos;
+ for (int i = 0; i < points.size()/3; i++)
+ {
+ double helper=0;
+ for (int j=0; j < points_neightri[i]->size();j++) helper += mybuf[m + points_neightri[i]->get(j)*size_one];
+ helper /= points_neightri[i]->size();
+ fprintf(fp,"%f\n",helper);
+ }
+ buf_pos++;
+ }
+ return;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void DumpMeshVTK::write_data_ascii_face(int n, double *mybuf)
{
int k, m, buf_pos;
diff --git a/src/dump_mesh_vtk.h b/src/dump_mesh_vtk.h
index dfe32625..3f50d4cc 100644
--- a/src/dump_mesh_vtk.h
+++ b/src/dump_mesh_vtk.h
@@ -44,6 +44,8 @@ class DumpMeshVTK : public Dump {
private: // column labels
+ int dataMode_;
+
int nMesh_;
class TriMesh **meshList_;
@@ -78,6 +80,8 @@ class DumpMeshVTK : public Dump {
void write_header_ascii(bigint ndump);
void write_data_ascii(int n, double *mybuf);
+ void write_data_ascii_point(int n, double *mybuf);
+ void write_data_ascii_face(int n, double *mybuf);
};
}
diff --git a/src/fix_ave_euler.cpp b/src/fix_ave_euler.cpp
new file mode 100644
index 00000000..65f86c01
--- /dev/null
+++ b/src/fix_ave_euler.cpp
@@ -0,0 +1,479 @@
+/* ----------------------------------------------------------------------
+ LIGGGHTS - LAMMPS Improved for General Granular and Granular Heat
+ Transfer Simulations
+
+ LIGGGHTS is part of the CFDEMproject
+ www.liggghts.com | www.cfdem.com
+
+ Christoph Kloss, christoph.kloss@cfdem.com
+ Copyright 2009-2012 JKU Linz
+ Copyright 2012- DCS Computing GmbH, Linz
+
+ LIGGGHTS is based on LAMMPS
+ LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+ http://lammps.sandia.gov, Sandia National Laboratories
+ Steve Plimpton, sjplimp@sandia.gov
+
+ This software is distributed under the GNU General Public License.
+
+ See the README file in the top-level directory.
+------------------------------------------------------------------------- */
+
+/* ----------------------------------------------------------------------
+ Contributing author for triclinic: Andreas Aigner (JKU)
+------------------------------------------------------------------------- */
+
+#include "stdlib.h"
+#include "string.h"
+#include "fix_ave_euler.h"
+#include "compute_stress_atom.h"
+#include "atom.h"
+#include "force.h"
+#include "domain.h"
+#include "modify.h"
+#include "neighbor.h"
+#include "update.h"
+#include "memory.h"
+#include "error.h"
+
+#define BIG 1000000000
+
+#define INVOKED_PERATOM 8
+
+using namespace LAMMPS_NS;
+using namespace FixConst;
+
+/* ---------------------------------------------------------------------- */
+
+FixAveEuler::FixAveEuler(LAMMPS *lmp, int narg, char **arg) :
+ Fix(lmp, narg, arg),
+ exec_every_(1),
+ box_change_(false),
+ cell_size_ideal_rel_(3.),
+ cell_size_ideal_(0.),
+ ncells_(0),
+ ncells_max_(0),
+ ncellptr_max_(0),
+ cellhead_(NULL),
+ cellptr_(NULL),
+ center_(NULL),
+ v_av_(NULL),
+ vol_fr_(NULL),
+ mass_(NULL),
+ stress_(NULL),
+ compute_stress_(NULL)
+{
+ // this fix produces a global array
+
+ array_flag = 1;
+ size_array_rows = BIG;
+ size_array_cols = 3 + 1 + 3;
+
+ triclinic_ = domain->triclinic; //NP modified A.A.
+
+ // parse args
+ if (narg < 6) error->all(FLERR,"Illegal fix ave/pic command");
+ int iarg = 3;
+
+ if(strcmp(arg[iarg++],"nevery"))
+ error->fix_error(FLERR,this,"expecting keyword 'nevery'");
+ exec_every_ = force->inumeric(arg[iarg++]);
+ if(exec_every_ < 1)
+ error->fix_error(FLERR,this,"'nevery' > 0 required");
+ nevery = exec_every_;
+
+ if(strcmp(arg[iarg++],"cell_size_relative"))
+ error->fix_error(FLERR,this,"expecting keyword 'cell_size_relative'");
+ cell_size_ideal_rel_ = force->numeric(arg[iarg++]);
+ if(cell_size_ideal_rel_ < 3.)
+ error->fix_error(FLERR,this,"'cell_size_relative' > 3 required");
+
+ //NP from FixPrint
+ // add nfirst to all computes that store invocation times
+ // since don't know a priori which are invoked via variables by this fix
+ // once in end_of_step() can set timestep for ones actually invoked
+
+ int nfirst = (update->ntimestep/nevery)*nevery + nevery;
+ modify->addstep_compute_all(nfirst);
+}
+
+/* ---------------------------------------------------------------------- */
+
+FixAveEuler::~FixAveEuler()
+{
+ memory->destroy(cellhead_);
+ memory->destroy(cellptr_);
+ memory->destroy(center_);
+ memory->destroy(v_av_);
+ memory->destroy(vol_fr_);
+ memory->destroy(mass_);
+ memory->destroy(stress_);
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixAveEuler::post_create()
+{
+ // register per particle properties
+ if(!compute_stress_)
+ {
+ char* arg[4];
+ arg[0]="stress_faveu";
+ arg[1]="all";
+ arg[2]="stress/atom";
+ arg[3]="pair";
+
+ modify->add_compute(4,arg);
+ compute_stress_ = static_cast(modify->compute[modify->find_compute(arg[0])]);
+ }
+}
+
+/* ---------------------------------------------------------------------- */
+
+int FixAveEuler::setmask()
+{
+ int mask = 0;
+ mask |= END_OF_STEP;
+ return mask;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixAveEuler::init()
+{
+ if(!atom->radius_flag)
+ error->fix_error(FLERR,this,"requires atom attribute radius");
+ if(!atom->rmass_flag)
+ error->fix_error(FLERR,this,"requires atom attribute mass");
+
+ // check if box constant
+ box_change_ = false;
+ if(domain->box_change)
+ box_change_ = true;
+}
+
+/* ----------------------------------------------------------------------
+ setup initial bins
+ only does averaging if nvalid = current timestep
+------------------------------------------------------------------------- */
+
+void FixAveEuler::setup(int vflag)
+{
+ setup_bins();
+ //end_of_step();
+}
+
+/* ----------------------------------------------------------------------
+ setup 3d bins and their extent and coordinates
+ called at setup() and when averaging occurs if box size changes
+ similar to FixAveSpatial::setup_bins() and PairDSMC::init_style()
+
+ bins are subbox - skin/2 so owned particles cannot move outside
+ bins - so do not have to extrapolate
+------------------------------------------------------------------------- */
+
+void FixAveEuler::setup_bins()
+{
+ int ibin;
+ double cell_size[3];
+
+ // calc ideal cell size as multiple of max cutoff
+ cell_size_ideal_ = cell_size_ideal_rel_ * (neighbor->cutneighmax-neighbor->skin);
+
+ cell_size[0] = cell_size[1] = cell_size[2] = cell_size_ideal_;
+
+ //NP TODO: There are lot if(triclinic_); Maybe a template function?
+ if (triclinic_) {
+ domain->x2lamda(cell_size,cell_size);
+ }
+
+ for(int dim = 0; dim < 3; dim++)
+ {
+ // get bounds
+ if (triclinic_ == 0) {
+ lo_[dim] = domain->sublo[dim];
+ hi_[dim] = domain->subhi[dim];
+ } else {
+ lo_lamda_[dim] = domain->sublo_lamda[dim];
+ hi_lamda_[dim] = domain->subhi_lamda[dim];
+ }
+ }
+ if (triclinic_) {
+ domain->lamda2x(lo_lamda_,lo_);
+ domain->lamda2x(hi_lamda_,hi_);
+ }
+
+ for(int dim = 0; dim < 3; dim++) {
+ // round down (makes cell size larger)
+ if (triclinic_) {
+ ncells_dim_[dim] = static_cast((hi_lamda_[dim]-lo_lamda_[dim])/cell_size[dim]);
+ cell_size_lamda_[dim] = (hi_lamda_[dim]-lo_lamda_[dim])/static_cast(ncells_dim_[dim]);
+ } else {
+ ncells_dim_[dim] = static_cast((hi_[dim]-lo_[dim])/cell_size_ideal_);
+ cell_size_[dim] = (hi_[dim]-lo_[dim])/static_cast(ncells_dim_[dim]);
+ }
+ }
+ if (triclinic_) domain->lamda2x(cell_size_lamda_,cell_size_);
+
+ for(int dim = 0; dim < 3; dim++)
+ {
+ cell_size_inv_[dim] = 1./cell_size_[dim];
+ if (triclinic_) cell_size_lamda_inv_[dim] = 1./cell_size_lamda_[dim];
+
+ // add 2 extra cells for ghosts
+ // do not need to match subbox expaned by cutneighmax
+ // cell size will always be larger than cutneighmax
+ //ncells_dim_[dim] += 2;
+
+ // do not change cell_size_ but adapt lo_/hi_
+ //lo_[dim] -= cell_size_[dim];
+ //hi_[dim] += cell_size_[dim];
+ }
+
+ ncells_ = ncells_dim_[0]*ncells_dim_[1]*ncells_dim_[2];
+ cell_volume_ = cell_size_[0]*cell_size_[1]*cell_size_[2];
+
+ // (re) allocate spatial bin arrays
+ if (ncells_ > ncells_max_)
+ {
+ ncells_max_ = ncells_;
+ memory->grow(cellhead_,ncells_max_,"ave/euler:cellhead_");
+ memory->grow(center_,ncells_max_,3,"ave/euler:center_");
+ memory->grow(v_av_, ncells_max_,3,"ave/euler:v_av_");
+ memory->grow(vol_fr_,ncells_max_, "ave/euler:vol_fr_");
+ memory->grow(mass_,ncells_max_, "ave/euler:mass_");
+ memory->grow(stress_,ncells_max_,7,"ave/euler:stress_");
+ }
+
+ // calculate center corrdinates for cells
+ for(int ix = 0; ix < ncells_dim_[0]; ix++)
+ {
+ for(int iy = 0; iy < ncells_dim_[1]; iy++)
+ {
+ for(int iz = 0; iz < ncells_dim_[2]; iz++)
+ {
+ ibin = iz*ncells_dim_[1]*ncells_dim_[0] + iy*ncells_dim_[0] + ix;
+
+ if (triclinic_) {
+ center_[ibin][0] = lo_lamda_[0] + (static_cast(ix)+0.5) * cell_size_lamda_[0];
+ center_[ibin][1] = lo_lamda_[1] + (static_cast(iy)+0.5) * cell_size_lamda_[1];
+ center_[ibin][2] = lo_lamda_[2] + (static_cast(iz)+0.5) * cell_size_lamda_[2];
+ domain->lamda2x(center_[ibin],center_[ibin]);
+
+ } else {
+ center_[ibin][0] = lo_[0] + (static_cast(ix)+0.5) * cell_size_[0];
+ center_[ibin][1] = lo_[1] + (static_cast(iy)+0.5) * cell_size_[1];
+ center_[ibin][2] = lo_[2] + (static_cast(iz)+0.5) * cell_size_[2];
+ }
+ }
+ }
+ }
+
+ // print to screen and log
+ /*
+ if (comm->me == 0)
+ {
+ if (screen) fprintf(screen,"Euler cell size on proc %d = %f (%d x %d x %d grid)\n",
+ comm->me,cell_size_ideal_,ncells_dim_[0],ncells_dim_[1],ncells_dim_[2]);
+ if (logfile) fprintf(logfile,"Euler cell size on proc %d = %f (%d x %d x %d grid)\n",
+ comm->me,cell_size_ideal_,ncells_dim_[0],ncells_dim_[1],ncells_dim_[2]);
+ }
+ */
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixAveEuler::end_of_step()
+{
+ /*NP do not need this for end_of_step()
+ if (update->ntimestep % exec_every_)
+ return;
+ */
+
+ // have to adapt grid if box changes
+ if(box_change_)
+ setup_bins();
+
+ // bin atoms
+ bin_atoms();
+
+ // calculate Eulerian grid properties
+ calculate_eu();
+}
+
+/* ----------------------------------------------------------------------
+ bin owned and ghost atoms
+ this also implies we do not need to wrap around PBCs
+ bin ghost atoms only if inside my grid
+------------------------------------------------------------------------- */
+
+void FixAveEuler::bin_atoms()
+{
+ int i,ibin;
+ double **x = atom->x;
+ int nall = atom->nlocal + atom->nghost;
+
+ for (i = 0; i < ncells_max_; i++)
+ cellhead_[i] = -1;
+
+ // re-alloc cellptr_ if necessary
+ if(nall > ncellptr_max_)
+ {
+ ncellptr_max_ = nall;
+ memory->grow(cellptr_,ncellptr_max_,"ave/pic:cellptr_");
+ }
+
+ // bin in reverse order so linked list will be in forward order
+ // also puts ghost atoms at end of list
+
+ for (i = nall-1; i >= 0; i--)
+ {
+ ibin = coord2bin(x[i]);
+
+ // ghosts outside grid may return -1
+ // lets ignore them
+ if (ibin < 0) {
+ //NP Test-output if any particle is out of bounds; modified A.A.
+ /*NL*/ //int nlocal = atom->nlocal;
+ /*NL*/ //if(screen) fprintf(screen,"Particle with id= %d is out of bounds. nlocal = %d\n",i,nlocal);
+ continue;
+ }
+
+ cellptr_[i] = cellhead_[ibin];
+ cellhead_[ibin] = i;
+ }
+}
+
+/* ----------------------------------------------------------------------
+ map coord to grid, also return ix,iy,iz indices in each dim
+------------------------------------------------------------------------- */
+
+inline int FixAveEuler::coord2bin(double *x)
+{
+ int i,iCell[3];
+ double float_iCell[3],tmp_x[3];
+
+ if (triclinic_) {
+ double tmp_x[3];
+ domain->x2lamda(x,tmp_x);
+ for (i=0;i<3;i++) {//NP TODO: Dimension?!
+ float_iCell[i] = (tmp_x[i]-lo_lamda_[i])*cell_size_lamda_inv_[i];
+ iCell[i] = static_cast (float_iCell[i] >= 0 ? float_iCell[i] : float_iCell[i]-1);
+ }
+ } else {
+ for (i=0;i<3;i++) {//NP TODO: Dimension?!
+ float_iCell[i] = (x[i]-lo_[i])*cell_size_inv_[i];
+ iCell[i] = static_cast (float_iCell[i] >= 0 ? float_iCell[i] : float_iCell[i]-1);
+ }
+ }
+
+ //NP Test-output if any particle is out of bounds; modified A.A.
+ /*NL*/ //if (iCell[0]>=ncells_dim_[0] || iCell[1]>=ncells_dim_[1] || iCell[2]>=ncells_dim_[2] || iCell[0]<0 || iCell[1]<0 || iCell[2]<0) {
+ /*NL*/ // if(screen) fprintf(screen,"Particle with x[0-2] = %f %f %f would have iCell[0] = %d, iCell[1] = %d and iCell[2] = %d \n",x[0],x[1],x[2],iCell[0],iCell[1],iCell[2]);
+ /*NL*/ // return -1;
+ /*NL*/ //}
+
+ return iCell[2]*ncells_dim_[1]*ncells_dim_[0] + iCell[1]*ncells_dim_[0] + iCell[0];
+}
+
+/* ----------------------------------------------------------------------
+ calculate Eulerian data, use interpolation function
+------------------------------------------------------------------------- */
+
+void FixAveEuler::calculate_eu()
+{
+ int ncount;
+ double **v = atom->v;
+ double *radius = atom->radius;
+ double *rmass = atom->rmass;
+
+ double prefactor_vol_fr = 4./3.*M_PI/cell_volume_;
+ double prefactor_stress = 1./cell_volume_;
+ double vel_x_mass[3];
+
+ // wrap compute with clear/add
+ modify->clearstep_compute();
+
+
+ /*NL*/ //fprintf(screen,"ts %d, exec_every_ %d\n",update->ntimestep,exec_every_);
+
+ // invoke compute if not previously invoked
+ //NP also mess around with invoked_flag b/c addstep_compute
+ //NP will only add a new step to those computes which were
+ //NP invoked this time-step
+ if (!(compute_stress_->invoked_flag & INVOKED_PERATOM)) {
+ compute_stress_->compute_peratom();
+ compute_stress_->invoked_flag |= INVOKED_PERATOM;
+ }
+
+ // need to get pointer here since compute_peratom() may realloc
+ double **stress_atom = compute_stress_->array_atom;
+
+ // loop all binned particles
+ // each particle can contribute to the cell that it has been binned
+ // to plus its 26 neighs
+
+ for(int icell = 0; icell < ncells_; icell++)
+ {
+ ncount = 0;
+ vectorZeroize3D(v_av_[icell]);
+ vol_fr_[icell] = 0.;
+ mass_[icell] = 0.;
+ vectorZeroizeN(stress_[icell],7);
+
+ // add contributions of particle - v and volume fraction
+ // v is favre-averaged (mass-averaged)
+
+ for(int iatom = cellhead_[icell/*+stencil*/]; iatom >= 0; iatom = cellptr_[iatom])
+ {
+ vectorScalarMult3D(v[iatom],rmass[iatom],vel_x_mass);
+ vectorAdd3D(v_av_[icell],vel_x_mass,v_av_[icell]);
+ vol_fr_[icell] += radius[iatom]*radius[iatom]*radius[iatom];
+ mass_[icell] += rmass[iatom];
+ ncount++;
+ }
+
+ /*NL*/ //if(ncount) fprintf(screen,"Cell %d has %d particles. \n",icell,ncount); // modified A.A.
+
+ if(ncount) vectorScalarDiv3D(v_av_[icell],static_cast(ncount)*mass_[icell]);
+ vol_fr_[icell] *= prefactor_vol_fr;
+
+ // add contribution of particle - stress
+ // need v before can calculate stress
+ // stress is molecular diffusion + contact forces
+
+ for(int iatom = cellhead_[icell/*+stencil*/]; iatom >= 0; iatom = cellptr_[iatom])
+ {
+ stress_[icell][1] += -rmass[iatom]*(v[iatom][0]-v_av_[icell][0])*(v[iatom][0]-v_av_[icell][0]) + stress_atom[iatom][0];
+ stress_[icell][2] += -rmass[iatom]*(v[iatom][1]-v_av_[icell][1])*(v[iatom][1]-v_av_[icell][1]) + stress_atom[iatom][1];
+ stress_[icell][3] += -rmass[iatom]*(v[iatom][2]-v_av_[icell][2])*(v[iatom][2]-v_av_[icell][2]) + stress_atom[iatom][2];
+ stress_[icell][4] += -rmass[iatom]*(v[iatom][0]-v_av_[icell][0])*(v[iatom][1]-v_av_[icell][1]) + stress_atom[iatom][3];
+ stress_[icell][5] += -rmass[iatom]*(v[iatom][0]-v_av_[icell][0])*(v[iatom][2]-v_av_[icell][2]) + stress_atom[iatom][4];
+ stress_[icell][6] += -rmass[iatom]*(v[iatom][1]-v_av_[icell][1])*(v[iatom][2]-v_av_[icell][2]) + stress_atom[iatom][5]; //NP modified A.A.: it was stress_[icell][5]
+ stress_[icell][0] = -0.333333333333333*(stress_[icell][1]+stress_[icell][2]+stress_[icell][3]);
+ }
+ vectorScalarMultN(7,stress_[icell],prefactor_stress);
+ }
+
+ // wrap with clear/add
+ modify->addstep_compute(update->ntimestep + exec_every_);
+}
+
+/* ----------------------------------------------------------------------
+ return I,J array value
+ if I exceeds current bins, return 0.0 instead of generating an error
+ column 1,2,3 = bin coords, next column = vol fr,
+ remaining columns = vel, stress
+------------------------------------------------------------------------- */
+
+double FixAveEuler::compute_array(int i, int j)
+{
+ if(i >= ncells_) return 0.0;
+
+ else if(j < 3) return center_[i][j];
+ else if(j == 3) return vol_fr_[i];
+ else if(j < 7) return v_av_[i][j-4];
+ else if(j == 7) return stress_[i][0];
+ else if(j < 14) return stress_[i][j-8];
+ else return 0.0;
+}
diff --git a/src/fix_ave_euler.h b/src/fix_ave_euler.h
new file mode 100644
index 00000000..80496e6a
--- /dev/null
+++ b/src/fix_ave_euler.h
@@ -0,0 +1,146 @@
+/* ----------------------------------------------------------------------
+ LIGGGHTS - LAMMPS Improved for General Granular and Granular Heat
+ Transfer Simulations
+
+ LIGGGHTS is part of the CFDEMproject
+ www.liggghts.com | www.cfdem.com
+
+ Christoph Kloss, christoph.kloss@cfdem.com
+ Copyright 2009-2012 JKU Linz
+ Copyright 2012- DCS Computing GmbH, Linz
+
+ LIGGGHTS is based on LAMMPS
+ LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+ http://lammps.sandia.gov, Sandia National Laboratories
+ Steve Plimpton, sjplimp@sandia.gov
+
+ This software is distributed under the GNU General Public License.
+
+ See the README file in the top-level directory.
+------------------------------------------------------------------------- */
+
+/* ----------------------------------------------------------------------
+ Contributing author for triclinic: Andreas Aigner (JKU)
+------------------------------------------------------------------------- */
+
+#ifdef FIX_CLASS
+
+FixStyle(ave/euler,FixAveEuler)
+FixStyle(ave/euler/stress,FixAveEuler)
+
+#else
+
+#ifndef LMP_FIX_AVE_EULER_H
+#define LMP_FIX_AVE_EULER_H
+
+#include "stdio.h"
+#include "fix.h"
+#include "vector_liggghts.h"
+
+namespace LAMMPS_NS {
+
+class FixAveEuler : public Fix {
+
+ public:
+
+ FixAveEuler(class LAMMPS *, int, char **);
+ ~FixAveEuler();
+
+ void post_create();
+ int setmask();
+ void init();
+ void setup(int vflag);
+
+ void end_of_step();
+
+ double compute_array(int i, int j);
+
+ // inline access functions for cell based values
+
+ inline int ncells()
+ { return ncells_; }
+
+ inline double cell_center(int i, int j)
+ { return center_[i][j]; }
+
+ inline double cell_v_av(int i, int j)
+ { return v_av_[i][j]; }
+
+ inline double cell_vol_fr(int i)
+ { return vol_fr_[i]; }
+
+ inline double cell_pressure(int i)
+ { return stress_[i][0]; }
+
+ inline double cell_stress(int i,int j)
+ { return stress_[i][j+1]; }
+
+ private:
+
+ void setup_bins();
+ void bin_atoms();
+ void calculate_eu();
+ inline int coord2bin(double *x); //NP modified A.A.
+
+ int exec_every_;
+ bool box_change_;
+ int triclinic_; //NP modified A.A.
+
+ // desired cell size over max particle diameter
+ double cell_size_ideal_rel_;
+
+ // desired cell size
+ double cell_size_ideal_;
+
+ // number of cells
+ int ncells_, ncells_dim_[3];
+
+ // extent of grid in xyz
+ double lo_[3],hi_[3];
+ double lo_lamda_[3],hi_lamda_[3]; //NP modified A.A.
+
+ // cell size and inverse size in xyz, cell and volume
+ double cell_size_[3];
+ double cell_size_inv_[3];
+ double cell_volume_;
+ double cell_size_lamda_[3]; //NP modified A.A.
+ double cell_size_lamda_inv_[3]; //NP modified A.A.
+
+ // length of cellhead_, center_, v_av_, vol_fr_ arrays
+ int ncells_max_;
+
+ // length of cellptr_ array
+ int ncellptr_max_;
+
+ // atom - cell mapping
+ int *cellhead_; // ptr to 1st atom in each cell
+ int *cellptr_; // ptr to next atom in each bin
+
+ /* --------- DATA --------- */
+
+ // cell center
+ double **center_;
+
+ // cell-based averaged velocity
+ double **v_av_;
+
+ // cell-based volume fraction
+ double *vol_fr_;
+
+ // cell-based mass
+ double *mass_;
+
+ // cell-based stress
+ // [0]: pressure
+ // [1-3]: 00-11-22
+ // [4-6]: 01-02-12
+ double **stress_;
+
+ // stress computation
+ class ComputeStressAtom *compute_stress_;
+};
+
+}
+
+#endif
+#endif
diff --git a/src/fix_insert.cpp b/src/fix_insert.cpp
index 09ed4265..21512aea 100644
--- a/src/fix_insert.cpp
+++ b/src/fix_insert.cpp
@@ -154,6 +154,13 @@ FixInsert::FixInsert(LAMMPS *lmp, int narg, char **arg) :
else error->fix_error(FLERR,this,"");
iarg += 2;
hasargs = true;
+ } else if (strcmp(arg[iarg],"verbose") == 0) {
+ if (iarg+2 > narg) error->fix_error(FLERR,this,"");
+ if(strcmp(arg[iarg+1],"no")==0) print_stats_during_flag = 0;
+ else if(strcmp(arg[iarg+1],"yes")==0) print_stats_during_flag = 1;
+ else error->fix_error(FLERR,this,"");
+ iarg += 2;
+ hasargs = true;
} else if (strcmp(arg[iarg],"vel") == 0) {
if (iarg+5 > narg) error->fix_error(FLERR,this,"");
if (strcmp(arg[iarg+1],"constant") == 0)
@@ -295,6 +302,8 @@ void FixInsert::init_defaults()
//NP initialize as unit maxtrix
quatUnitize4D(quat_insert);
+
+ print_stats_during_flag = 1;
}
/* ---------------------------------------------------------------------- */
@@ -316,11 +325,11 @@ void FixInsert::print_stats_start()
{
if (me == 0 && print_stats_start_flag && ninsert_exists) {
if (screen)
- fprintf(screen ,"Particle insertion: %f particles every %d steps - particle rate %f (mass rate %f)\n %d particles (mass %f) within %d steps\n",
+ fprintf(screen ,"INFO: Particle insertion: %f particles every %d steps - particle rate %f (mass rate %f)\n %d particles (mass %f) within %d steps\n",
ninsert_per,insert_every,nflowrate,massflowrate,ninsert,massinsert,final_ins_step-first_ins_step);
if (logfile)
- fprintf(logfile,"Particle insertion: %f particles every %d steps - particle rate %f, (mass rate %f)\n %d particles (mass %f) within %d steps\n",
+ fprintf(logfile,"INFO: Particle insertion: %f particles every %d steps - particle rate %f, (mass rate %f)\n %d particles (mass %f) within %d steps\n",
ninsert_per,insert_every,nflowrate,massflowrate,ninsert,massinsert,final_ins_step-first_ins_step);
}
}
@@ -331,14 +340,14 @@ void FixInsert::print_stats_during(int ninsert_this, double mass_inserted_this)
{
int step = update->ntimestep;
- if (me == 0)
+ if (me == 0 && print_stats_during_flag)
{
if (screen)
- fprintf(screen ,"Particle insertion: inserted %d particle templates (mass %f) at step %d\n - a total of %d particle templates (mass %f) inserted so far.\n",
+ fprintf(screen ,"INFO: Particle insertion: inserted %d particle templates (mass %f) at step %d\n - a total of %d particle templates (mass %f) inserted so far.\n",
ninsert_this,mass_inserted_this,step,ninserted,massinserted);
if (logfile)
- fprintf(logfile,"Particle insertion: inserted %d particle templates (mass %f) at step %d\n - a total of %d particle templates (mass %f) inserted so far.\n",
+ fprintf(logfile,"INFO: Particle insertion: inserted %d particle templates (mass %f) at step %d\n - a total of %d particle templates (mass %f) inserted so far.\n",
ninsert_this,mass_inserted_this,step,ninserted,massinserted);
}
}
diff --git a/src/fix_insert.h b/src/fix_insert.h
index 4e4791d8..d7be7bd8 100644
--- a/src/fix_insert.h
+++ b/src/fix_insert.h
@@ -136,6 +136,7 @@ class FixInsert : public Fix {
// determine if print stats
int print_stats_start_flag;
+ int print_stats_during_flag;
class FixMultisphere *fix_multisphere;
class Multisphere *multisphere;
diff --git a/src/fix_insert_rate_region.cpp b/src/fix_insert_rate_region.cpp
index 2020d469..a8492b0c 100644
--- a/src/fix_insert_rate_region.cpp
+++ b/src/fix_insert_rate_region.cpp
@@ -51,6 +51,15 @@ FixInsertRateRegion::FixInsertRateRegion(LAMMPS *lmp, int narg, char **arg) :
{
// no fixed total number of particles inserted by this fix exists
ninsert_exists = 1;
+
+ bool hasargs = true;
+ while(iarg < narg && hasargs)
+ {
+ hasargs = false;
+ if (strcmp(arg[iarg],"some_arg") == 0) {
+ } else if(strcmp(style,"insert/rate/region") == 0)
+ error->fix_error(FLERR,this,"unknown keyword");
+ }
}
/* ---------------------------------------------------------------------- */
diff --git a/src/fix_mesh.h b/src/fix_mesh.h
index 1f6aed50..8c566634 100644
--- a/src/fix_mesh.h
+++ b/src/fix_mesh.h
@@ -36,7 +36,6 @@
#include "fix.h"
-
namespace LAMMPS_NS
{
class FixMesh : public Fix
diff --git a/src/fix_mesh_surface.cpp b/src/fix_mesh_surface.cpp
index 6d9fc96e..79db4cb8 100644
--- a/src/fix_mesh_surface.cpp
+++ b/src/fix_mesh_surface.cpp
@@ -30,6 +30,7 @@
#include
#include
#include "error.h"
+#include "group.h"
#include "force.h"
#include "bounding_box.h"
#include "input_mesh_tri.h"
@@ -197,7 +198,7 @@ void FixMeshSurface::final_integrate()
called from fix wall/gran out of post_create()
------------------------------------------------------------------------- */
-void FixMeshSurface::createNeighList()
+void FixMeshSurface::createNeighList(int igrp)
{
if(fix_mesh_neighlist_) return;
char *neighlist_name = new char[strlen(id)+1+10];
@@ -213,6 +214,10 @@ void FixMeshSurface::createNeighList()
fix_mesh_neighlist_ =
static_cast(modify->find_fix_id(neighlist_name));
+ // fix added with "all", correct this now
+ fix_mesh_neighlist_->igroup = igrp;
+ fix_mesh_neighlist_->groupbit = group->bitmask[igrp];
+
delete []fixarg;
delete []neighlist_name;
diff --git a/src/fix_mesh_surface.h b/src/fix_mesh_surface.h
index 4d921673..2c725964 100644
--- a/src/fix_mesh_surface.h
+++ b/src/fix_mesh_surface.h
@@ -60,7 +60,7 @@ namespace LAMMPS_NS
virtual void pre_force(int);
virtual void final_integrate();
- void createNeighList();
+ void createNeighList(int igrp);
void createContactHistory(int dnum);
inline bool trackStress()
diff --git a/src/fix_neighlist_mesh.cpp b/src/fix_neighlist_mesh.cpp
index 6a6d094f..313e5bcb 100644
--- a/src/fix_neighlist_mesh.cpp
+++ b/src/fix_neighlist_mesh.cpp
@@ -190,6 +190,8 @@ void FixNeighlistMesh::pre_force(int vflag)
void FixNeighlistMesh::handleTriangle(int iTri)
{
+ int *mask = atom->mask;
+
// get bounding box of element on this subdomain
//NP is bounded by sublo, subhi
BoundingBox b = mesh_->getElementBoundingBoxOnSubdomain(iTri);
@@ -235,6 +237,8 @@ void FixNeighlistMesh::handleTriangle(int iTri)
int iAtom = binhead[iBin];
while(iAtom != -1 && iAtom < nlocal)
{
+ if(! (mask[iAtom] & groupbit)) continue;
+
/*NL*/ if(DEBUGMODE_LMP_FIX_NEIGHLIST_MESH && DEBUG_LMP_FIX_NEIGHLIST_MESH_M_ID == mesh_->id(iTri)
/*NL*/ && DEBUG_LMP_FIX_NEIGHLIST_MESH_P_ID == atom->map(iAtom) )
/*NL*/ fprintf(screen, " handleTriangle atom %d for tri id %d on proc %d\n",
diff --git a/src/fix_wall_gran.cpp b/src/fix_wall_gran.cpp
index bf031727..b0d13aae 100644
--- a/src/fix_wall_gran.cpp
+++ b/src/fix_wall_gran.cpp
@@ -284,7 +284,7 @@ void FixWallGran::post_create()
// also create contact tracker
for(int i=0;icreateNeighList();
+ FixMesh_list_[i]->createNeighList(igroup);
if(dnum()>0)FixMesh_list_[i]->createContactHistory(dnum());
}
}
@@ -485,6 +485,8 @@ void FixWallGran::post_force_respa(int vflag, int ilevel, int iloop)
void FixWallGran::post_force_mesh(int vflag)
{
+ //NP groupbit accounted for via neighlist/mesh
+
// contact properties
double force_old[3],force_wall[3],v_wall[3],bary[3];
double delta[3],deltan;
@@ -551,8 +553,8 @@ void FixWallGran::post_force_mesh(int vflag)
/*NL*/// if(!radius_) error->one(FLERR,"need to revise this for SPH");
deltan = mesh->resolveTriSphereContactBary(iTri,radius_ ? radius_[iPart]:r0_ ,x_[iPart],delta,bary);
- /*NL*/// if(comm->me == 3 && update->ntimestep == 3735)// && mesh->id(iTri) == 4942)
- /*NL*/// fprintf(screen,"proc 3 moving mesh %d Tri %d iCont %d numNeigh %d 2, deltan %f\n",iMesh,mesh->id(iTri),iCont,numNeigh[iTri],deltan);
+ /*NL*/// if(atom->tag[iPart] == 624)// && mesh->id(iTri) == 4942)
+ /*NL*/// fprintf(screen,"particle 624 step %d moving mesh %d Tri %d iCont %d numNeigh %d 2, deltan %f\n",update->ntimestep,iMesh,mesh->id(iTri),iCont,numNeigh[iTri],deltan);
if(deltan > 0.)
{
@@ -569,7 +571,7 @@ void FixWallGran::post_force_mesh(int vflag)
for(int i = 0; i < 3; i++)
v_wall[i] = (bary[0]*vMesh[iTri][0][i] + bary[1]*vMesh[iTri][1][i] + bary[2]*vMesh[iTri][2][i]);
- post_force_eval_contact(iPart,deltan,delta,v_wall,c_history,FixMesh_list_[iMesh],mesh,iTri);
+ post_force_eval_contact(iPart,deltan,delta,v_wall,c_history,iMesh,FixMesh_list_[iMesh],mesh,iTri);
}
/*NL*/// if(comm->me == 3 && update->ntimestep == 3735)// && mesh->id(iTri) == 4942)
@@ -600,6 +602,9 @@ void FixWallGran::post_force_mesh(int vflag)
int idTri = mesh->id(iTri);
deltan = mesh->resolveTriSphereContact(iTri,radius_ ? radius_[iPart]:r0_,x_[iPart],delta);
+ /*NL*/// if(atom->tag[iPart] == 624)// && mesh->id(iTri) == 4942)
+ /*NL*/// fprintf(screen,"particle 624 step %d nonmoving mesh %d Tri %d iCont %d numNeigh %d 2, deltan %f\n",update->ntimestep,iMesh,mesh->id(iTri),iCont,numNeigh[iTri],deltan);
+
/*NL*/// if(DEBUGMODE_LMP_FIX_WALL_GRAN && DEBUG_LMP_FIX_FIX_WALL_GRAN_M_ID == mesh->id(iTri) &&
/*NL*/// DEBUG_LMP_FIX_FIX_WALL_GRAN_P_ID == atom->tag[iPart])
/*NL*/// fprintf(screen,"step %d: handling tri id %d with particle id %d, deltan %f\n",update->ntimestep,mesh->id(iTri),atom->tag[iPart],deltan);
@@ -612,7 +617,7 @@ void FixWallGran::post_force_mesh(int vflag)
else //NP always handle contact for SPH
{
if(fix_contact && ! fix_contact->handleContact(iPart,idTri,c_history)) continue;
- post_force_eval_contact(iPart,deltan,delta,v_wall,c_history,FixMesh_list_[iMesh],mesh,iTri);
+ post_force_eval_contact(iPart,deltan,delta,v_wall,c_history,iMesh,FixMesh_list_[iMesh],mesh,iTri);
}
}
}
@@ -633,6 +638,8 @@ void FixWallGran::post_force_mesh(int vflag)
void FixWallGran::post_force_primitive(int vflag)
{
+ int *mask = atom->mask;
+
// contact properties
double force_old[3],force_wall[3];
double delta[3],deltan;
@@ -649,6 +656,9 @@ void FixWallGran::post_force_primitive(int vflag)
for (int iCont = 0; iCont < nNeigh ; iCont++, neighborList++)
{
int iPart = *neighborList;
+
+ if(!(mask[iPart] & groupbit)) continue;
+
deltan = primitiveWall_->resolveContact(x_[iPart],radius_?radius_[iPart]:r0_,delta);
/*NL*/ if(DEBUGMODE_LMP_FIX_WALL_GRAN && DEBUG_LMP_FIX_FIX_WALL_GRAN_P_ID == atom->tag[iPart])
@@ -671,7 +681,7 @@ void FixWallGran::post_force_primitive(int vflag)
------------------------------------------------------------------------- */
inline void FixWallGran::post_force_eval_contact(int iPart, double deltan, double *delta,
- double *v_wall, double *c_history, FixMeshSurface *fix_mesh, TriMesh *mesh, int iTri)
+ double *v_wall, double *c_history, int iMesh, FixMeshSurface *fix_mesh, TriMesh *mesh, int iTri)
{
double delr = (radius_ ? radius_[iPart] : r0_) + deltan;
@@ -714,6 +724,14 @@ inline void FixWallGran::post_force_eval_contact(int iPart, double deltan, doubl
}
}
+ // add to cwl
+ if(cwl_ && !addflag_)
+ {
+ double contactPoint[3];
+ vectorAdd3D(x_[iPart],delta,contactPoint);
+ cwl_->add_wall_1(iMesh,mesh->id(iTri),iPart,contactPoint);
+ }
+
// add heat flux
if(heattransfer_flag_)
addHeatFlux(mesh,iPart,deltan,1.);
diff --git a/src/fix_wall_gran.h b/src/fix_wall_gran.h
index c2679247..f4eea99d 100644
--- a/src/fix_wall_gran.h
+++ b/src/fix_wall_gran.h
@@ -158,7 +158,7 @@ class FixWallGran : public Fix {
void post_force_primitive(int);
inline void post_force_eval_contact(int iPart, double deltan, double *delta, double *v_wall,
- double *c_history, FixMeshSurface *fix_mesh = 0, TriMesh *mesh = 0, int iTri = 0);
+ double *c_history, int iMesh = 0, FixMeshSurface *fix_mesh = 0, TriMesh *mesh = 0, int iTri = 0);
};
}
diff --git a/src/multi_node_mesh_parallel_buffer_I.h b/src/multi_node_mesh_parallel_buffer_I.h
index e8b672e5..396e5b65 100644
--- a/src/multi_node_mesh_parallel_buffer_I.h
+++ b/src/multi_node_mesh_parallel_buffer_I.h
@@ -132,7 +132,7 @@
double *bufMesh = NULL, *sendbufElems = NULL, *recvbufElems = NULL;
bool dummy = false;
- /*NL*/ fprintf(this->screen,"buffersize mesh %d buffersize elems local %d\n",sizeMesh,sizeElements);
+ /*NL*/ //fprintf(this->screen,"buffersize mesh %d buffersize elems local %d\n",sizeMesh,sizeElements);
// pack global data into buffer
// do this only on proc 0
@@ -157,7 +157,7 @@
//NP send from all to proc 0
sizeElements_all = MPI_Gather0_Vector(sendbufElems,sizeElements,recvbufElems,this->world);
- /*NL*/ fprintf(this->screen,"buffersize mesh %d buffersize elems local %d, buffersize elems global %d\n",sizeMesh,sizeElements, sizeElements_all);
+ /*NL*/ //fprintf(this->screen,"buffersize mesh %d buffersize elems local %d, buffersize elems global %d\n",sizeMesh,sizeElements, sizeElements_all);
// actually write data to restart file
// do this only on proc 0
diff --git a/src/style_compute.h b/src/style_compute.h
index 9e129699..c7a0e31a 100644
--- a/src/style_compute.h
+++ b/src/style_compute.h
@@ -6,6 +6,7 @@
#include "compute_cna_atom.h"
#include "compute_com.h"
#include "compute_com_molecule.h"
+#include "compute_contact_atom.h"
#include "compute_coord_atom.h"
#include "compute_dihedral_local.h"
#include "compute_displace_atom.h"
diff --git a/src/style_dump.h b/src/style_dump.h
index e05bb602..74461452 100644
--- a/src/style_dump.h
+++ b/src/style_dump.h
@@ -3,6 +3,7 @@
#include "dump_custom.h"
#include "dump_dcd.h"
#include "dump_decomposition_vtk.h"
+#include "dump_euler_vtk.h"
#include "dump_image.h"
#include "dump_local.h"
#include "dump_mesh_stl.h"
diff --git a/src/style_fix.h b/src/style_fix.h
index 9207d9ef..90d977c8 100644
--- a/src/style_fix.h
+++ b/src/style_fix.h
@@ -2,6 +2,7 @@
#include "fix_addforce.h"
#include "fix_ave_atom.h"
#include "fix_ave_correlate.h"
+#include "fix_ave_euler.h"
#include "fix_ave_histo.h"
#include "fix_ave_spatial.h"
#include "fix_ave_time.h"
diff --git a/src/tri_mesh.h b/src/tri_mesh.h
index d29447e1..5c2015fa 100644
--- a/src/tri_mesh.h
+++ b/src/tri_mesh.h
@@ -48,6 +48,7 @@ namespace LAMMPS_NS
double resolveTriSphereContact(int nTri, double rSphere, double *cSphere, double *delta);
double resolveTriSphereContactBary(int nTri, double rSphere, double *cSphere, double *contactPoint,double *bary);
+ double resolveTriSphereContactBaryDefunct(int nTri, double rSphere, double *cSphere, double *contactPoint,double *bary);
bool resolveTriSphereNeighbuild(int nTri, double rSphere, double *cSphere, double treshold);
int generateRandomOwnedGhost(double *pos);
diff --git a/src/tri_mesh_I.h b/src/tri_mesh_I.h
index b189737c..a9034e04 100644
--- a/src/tri_mesh_I.h
+++ b/src/tri_mesh_I.h
@@ -30,8 +30,8 @@
#define SMALL_TRIMESH 1.e-10
-/*NL*/ #define DEBUGMODE_LMP_TRI_MESH_I_H false
-/*NL*/ #define DEBUGMODE_LMP_TRI_MESH_I_H_MESH_ID 21
+/*NL*/ #define DEBUGMODE_LMP_TRI_MESH_I_H false // (update->ntimestep > 3700)
+/*NL*/ #define DEBUGMODE_LMP_TRI_MESH_I_H_MESH_ID 5910
/* ---------------------------------------------------------------------- */
@@ -41,7 +41,7 @@
// coded in resolveTriSphereNeighbuild
/*NL*/ if(DEBUGMODE_LMP_TRI_MESH_I_H && DEBUGMODE_LMP_TRI_MESH_I_H_MESH_ID == id(nTri))
- /*NL*/ fprintf(screen, "step %d: resolveTriSphereContactNeigh for tri id %d with center %f %f %f\n",
+ /*NL*/ fprintf(screen, "step %d: resolveTriSphereContact for tri id %d with center %f %f %f\n",
/*NL*/ update->ntimestep,id(nTri),center_(nTri)[0],center_(nTri)[1],center_(nTri)[2]);
double tmp[3];
@@ -69,7 +69,8 @@
**edgeNorm = SurfaceMesh<3>::edgeNorm(nTri);
int i;
double distFromEdge(0.);
- for(i=0;i<3;i++){
+ for(i = 0; i < 3; i++)
+ {
vectorSubtract3D(csPlane,node[i],nodeToCsPlane);
distFromEdge = vectorDot3D(edgeNorm[i],nodeToCsPlane);
/*NL*/if(DEBUGMODE_LMP_TRI_MESH_I_H && DEBUGMODE_LMP_TRI_MESH_I_H_MESH_ID == id(nTri)) fprintf(screen,"distFromEdge %g\n",distFromEdge);
@@ -84,28 +85,59 @@
return (calcDist(cSphere,csPlane,delta) - rSphere);
}
- double *edgeVec = SurfaceMesh<3>::edgeVec(nTri)[i];
- double distFromNode = vectorDot3D(nodeToCsPlane,edgeVec);
+ double distFromNode;
- if(distFromNode < 0.)
+ // check for dist to next edge
+ // have to do this for obtuse angled triangle
+ bool isEdgeContNext = false;
+ if(i != 2)
{
- if(SurfaceMesh<3>::cornerActive(nTri)[i])
- {
- /*NL*/ if(DEBUGMODE_LMP_TRI_MESH_I_H && DEBUGMODE_LMP_TRI_MESH_I_H_MESH_ID == id(nTri)) fprintf(screen,"corner contact detected\n");
- return calcDist(cSphere,node[i],delta) - rSphere;
- }
- else
- return 1.;
+ int iNext = i+1;
+ double nodeToCsPlaneNext[3], distFromEdgeNext;
+ vectorSubtract3D(csPlane,node[iNext],nodeToCsPlaneNext);
+ distFromEdgeNext = vectorDot3D(edgeNorm[iNext],nodeToCsPlaneNext);
+ if(distFromEdgeNext > 0.)
+ {
+ double *edgeVecNext = SurfaceMesh<3>::edgeVec(nTri)[i+1];
+ distFromNode = vectorDot3D(nodeToCsPlaneNext,edgeVecNext);
+ if(distFromNode > -SMALL_TRIMESH)
+ {
+ isEdgeContNext = true;
+ i = iNext;
+ }
+ }
}
- else if(distFromNode > edgeLen(nTri)[i])
+
+ double *edgeVec = SurfaceMesh<3>::edgeVec(nTri)[i];
+
+ // go for corner contact only if not edge contact with next edge (case obtuse angled triangle)
+ if(!isEdgeContNext)
{
- if(SurfaceMesh<3>::cornerActive(nTri)[(i+1)%3])
- {
- /*NL*/ if(DEBUGMODE_LMP_TRI_MESH_I_H && DEBUGMODE_LMP_TRI_MESH_I_H_MESH_ID == id(nTri)) fprintf(screen,"corner contact detected\n");
- return calcDist(cSphere,node[(i+1)%3],delta) - rSphere;
- }
- else
- return 1.;
+ distFromNode = vectorDot3D(nodeToCsPlane,edgeVec);
+ /*NL*/ if(DEBUGMODE_LMP_TRI_MESH_I_H && DEBUGMODE_LMP_TRI_MESH_I_H_MESH_ID == id(nTri)) fprintf(screen,"distFromNode %g edgeLen(nTri)[i] %g\n",distFromNode,edgeLen(nTri)[i]);
+
+ if(distFromNode < 0.)
+ {
+ /*NL*/ if(DEBUGMODE_LMP_TRI_MESH_I_H && DEBUGMODE_LMP_TRI_MESH_I_H_MESH_ID == id(nTri)) fprintf(screen,"corner contact detected (0), but pot deactivated\n");
+ if(SurfaceMesh<3>::cornerActive(nTri)[i])
+ {
+ /*NL*/ if(DEBUGMODE_LMP_TRI_MESH_I_H && DEBUGMODE_LMP_TRI_MESH_I_H_MESH_ID == id(nTri)) fprintf(screen,"corner contact detected\n");
+ return calcDist(cSphere,node[i],delta) - rSphere;
+ }
+ else
+ return 1.;
+ }
+ else if(distFromNode > edgeLen(nTri)[i])
+ {
+ /*NL*/ if(DEBUGMODE_LMP_TRI_MESH_I_H && DEBUGMODE_LMP_TRI_MESH_I_H_MESH_ID == id(nTri)) fprintf(screen,"corner contact detected (1), but pot deactivated\n");
+ if(SurfaceMesh<3>::cornerActive(nTri)[(i+1)%3])
+ {
+ /*NL*/ if(DEBUGMODE_LMP_TRI_MESH_I_H && DEBUGMODE_LMP_TRI_MESH_I_H_MESH_ID == id(nTri)) fprintf(screen,"corner contact detected\n");
+ return calcDist(cSphere,node[(i+1)%3],delta) - rSphere;
+ }
+ else
+ return 1.;
+ }
}
/*NL*/ if(DEBUGMODE_LMP_TRI_MESH_I_H && DEBUGMODE_LMP_TRI_MESH_I_H_MESH_ID == id(nTri))
@@ -124,10 +156,41 @@
return d - rSphere;
}
+
/* ---------------------------------------------------------------------- */
inline double TriMesh::resolveTriSphereContactBary(int nTri, double rSphere,
double *cSphere, double *delta, double *bary)
+ {
+ double deltan, contactPoint[3], node0ToCsPlane[3];
+
+ // for now, use standard algorithm for distance calculation
+ deltan = resolveTriSphereContact(nTri, rSphere, cSphere, delta);
+
+ // additionaly calc bary coordinates
+ if(deltan < 0.)
+ {
+ vectorAdd3D(cSphere,delta,contactPoint);
+
+ double node0ToSphereCenter[3];
+ double *surfNorm = SurfaceMesh<3>::surfaceNorm(nTri);
+ vectorSubtract3D(cSphere,node_(nTri)[0],node0ToSphereCenter);
+
+ double csPlane[3],tmp[3];
+ double dNorm = vectorDot3D(surfNorm,node0ToSphereCenter);
+ vectorScalarMult3D(surfNorm,dNorm,tmp);
+ vectorSubtract3D(cSphere,tmp,csPlane);
+
+ MathExtraLiggghts::calcBaryTriCoords(node0ToCsPlane,edgeVec(nTri),edgeLen(nTri),bary);
+ }
+
+ return deltan;
+ }
+
+ /* ---------------------------------------------------------------------- */
+
+ inline double TriMesh::resolveTriSphereContactBaryDefunct(int nTri, double rSphere,
+ double *cSphere, double *delta, double *bary)
{
bool print = false;//(nTri == 25);
// this is only the overlap algorithm, neighbor list build is
@@ -142,9 +205,11 @@
double **n = node_(nTri);
// sphere-plane distance is coded explicitly because we need an intermediate result
+
double node0ToSphereCenter[3];
double *surfNorm = SurfaceMesh<3>::surfaceNorm(nTri);
vectorSubtract3D(cSphere,n[0],node0ToSphereCenter);
+
// normal distance of sphere center_ to plane
double dNorm = vectorDot3D(surfNorm,node0ToSphereCenter);
@@ -163,42 +228,35 @@
//NP optimized decision making
//NP use SMALL_TRIMESH to artificially enlarge triangle
int barySign = (bary[0] > -SMALL_TRIMESH) + 2*(bary[1] > -SMALL_TRIMESH) + 4*(bary[2] > -SMALL_TRIMESH);
-/*
- if(print){
- printf("node_ ");
- for(int i=0;i<3;i++)
- printf("%f %f %f | ",n[i][0],n[i][1],n[i][2]);
- printf("%f %f %f\n",csPlane[0],csPlane[1],csPlane[2]);
-
- printf("iTri %d barySign %d bary %f %f %f \n",nTri,barySign,bary[0],bary[1],bary[2]);
- }
-*/
double d(0.);
- switch(barySign){
- case 1:
- case 2:
- case 3: // bary[2] < 0 --> edge contact on edge[0]
- d = resolveEdgeContact(nTri,0,cSphere,csPlane,delta,bary);
- break;
- case 4:
- case 6: // bary[0] < 0 --> edge contact on edge[1]
- d = resolveEdgeContact(nTri,1,cSphere,csPlane,delta,bary);
- break;
- case 5: // bary[1] < 0 --> edge contact on edge[2]
- d = resolveEdgeContact(nTri,2,cSphere,csPlane,delta,bary);
- break;
- case 7: // face contact - all three barycentric coordinates are > 0
- d = calcDist(cSphere,csPlane,delta);
- break;
- default:
- d = 1.; // doesn't exist, just to satisfy the compiler
- break;
+ /*NL*/ if(DEBUGMODE_LMP_TRI_MESH_I_H && DEBUGMODE_LMP_TRI_MESH_I_H_MESH_ID == id(nTri))
+ /*NL*/ fprintf(screen,"step %d: detected barysign %d bary %f %f %f \n",update->ntimestep,barySign,bary[0],bary[1],bary[2]);
+
+ switch(barySign)
+ {
+ case 1:
+ case 2:
+ case 3: // bary[2] < 0 --> edge contact on edge[0]
+ d = resolveEdgeContact(nTri,0,cSphere,csPlane,delta,bary);
+ break;
+ case 4:
+ case 6: // bary[0] < 0 --> edge contact on edge[1]
+ d = resolveEdgeContact(nTri,1,cSphere,csPlane,delta,bary);
+ break;
+ case 5: // bary[1] < 0 --> edge contact on edge[2]
+ d = resolveEdgeContact(nTri,2,cSphere,csPlane,delta,bary);
+ break;
+ case 7: // face contact - all three barycentric coordinates are > 0
+ d = calcDist(cSphere,csPlane,delta);
+ break;
+ default:
+ d = 1.; // doesn't exist, just to satisfy the compiler
+ break;
}
return d == 1. ? d : d - rSphere;
-
}
/* ---------------------------------------------------------------------- */
@@ -217,17 +275,27 @@
double d(0.), distFromNode = vectorDot3D(tmp,edgeVec(iTri)[iEdge]);
if(distFromNode <= 0){
- if(!cornerActive(iTri)[iEdge]) d=1.;
- else{
- bary[iEdge] = 1.; bary[ip] = 0.; bary[ipp] = 0.;
- d = calcDist(p,node_(iTri)[iEdge],delta);
- }
+ double distFromPrevNode = vectorDot3D(pPlane,edgeVec(iTri)[ipp]);
+ if(distFromPrevNode < edgeLen(iTri)[ipp])
+ d = resolveEdgeContact(iTri,ipp,p,pPlane,delta,bary);
+ else{
+ if(!cornerActive(iTri)[iEdge]) d=1.;
+ else{
+ bary[iEdge] = 1.; bary[ip] = 0.; bary[ipp] = 0.;
+ d = calcDist(p,node_(iTri)[iEdge],delta);
+ }
+ }
} else if(distFromNode >= edgeLen(iTri)[iEdge]){
- if(!cornerActive(iTri)[ip]) d=1.;
- else{
- bary[iEdge] = 0.; bary[ip] = 1.; bary[ipp] = 0.;
- d = calcDist(p,node_(iTri)[ip],delta);
- }
+ double distFromNextNode = vectorDot3D(pPlane,edgeVec(iTri)[ip]);
+ if(distFromNextNode > 0)
+ d = resolveEdgeContact(iTri,ip,p,pPlane,delta,bary);
+ else{
+ if(!cornerActive(iTri)[ip]) d=1.;
+ else{
+ bary[iEdge] = 0.; bary[ip] = 1.; bary[ipp] = 0.;
+ d = calcDist(p,node_(iTri)[ip],delta);
+ }
+ }
} else{
if(!edgeActive(iTri)[iEdge]) d=1.;
else{
diff --git a/src/vector_liggghts.h b/src/vector_liggghts.h
index 02dcc2d4..ec464411 100644
--- a/src/vector_liggghts.h
+++ b/src/vector_liggghts.h
@@ -196,6 +196,12 @@ inline void vectorCopy4D(double *from,double *to)
to[3]=from[3];
}
+inline void vectorScalarMultN(int n,double *v,double s)
+{
+ for(int i = 0; i < n; i++)
+ v[i] = s*v[i];
+}
+
inline void vectorScalarMult3D(double *v,double s)
{
v[0]=s*v[0];
diff --git a/src/version_liggghts.h b/src/version_liggghts.h
index f9511bea..77b6ad5b 100644
--- a/src/version_liggghts.h
+++ b/src/version_liggghts.h
@@ -1 +1 @@
-#define LIGGGHTS_VERSION "LIGGGHTS-PFM 2.0.6, compiled 2012-08-08-11:10:32 by ckloss"
+#define LIGGGHTS_VERSION "LIGGGHTS-PFM 2.0.7, compiled 2012-08-23-17:21:27 by ckloss"
diff --git a/src/version_liggghts.txt b/src/version_liggghts.txt
index 157e54f3..f1547e6d 100644
--- a/src/version_liggghts.txt
+++ b/src/version_liggghts.txt
@@ -1 +1 @@
-2.0.6
+2.0.7