diff --git a/doc/Eqs/model_C_part_1.png b/doc/Eqs/model_C_part_1.png
index 6a18bdedafe93a7bcb68cbc7e64a4bf77a0e7280..ec6e131607758f8c8bbabe4465374016e01a4dd7 100644
Binary files a/doc/Eqs/model_C_part_1.png and b/doc/Eqs/model_C_part_1.png differ
diff --git a/doc/Eqs/model_C_part_2.png b/doc/Eqs/model_C_part_2.png
index 2bd4c4b89db12a24f911968a073a739b81c311c5..e8030e83f9fc05a514617347936f32051dfdcf71 100644
Binary files a/doc/Eqs/model_C_part_2.png and b/doc/Eqs/model_C_part_2.png differ
diff --git a/doc/Eqs/model_b1.png b/doc/Eqs/model_b1.png
index a44f9ff79a4acac9780ab09a5b03651b8d04a523..c993913ad10fc602723900613ab82233ce1b83c3 100644
Binary files a/doc/Eqs/model_b1.png and b/doc/Eqs/model_b1.png differ
diff --git a/doc/Manual.pdf b/doc/Manual.pdf
index 6fe72c9bb17c8efa4b2b77e48d4e9005a812f586..f86b7225de2c0b8471fde55790cc09bf00eb8a6c 100644
Binary files a/doc/Manual.pdf and b/doc/Manual.pdf differ
diff --git a/doc/set.html b/doc/set.html
index 6cfebeb9d440baa108b7721955cfa2b11edd9799..34a9b4964319e0dffaa747639ae8bed1f1ac8e80 100644
--- a/doc/set.html
+++ b/doc/set.html
@@ -68,6 +68,10 @@
<I>meso_e</I> value = energy of SPH particles (need units)
<I>meso_cv</I> value = heat capacity of SPH particles (need units)
<I>meso_rho</I> value = density of SPH particles (need units)
+ <I>add</I> value = yes no
+ yes = add per-atom quantities to a region or a group
+ <I>until</I> value = final timestep
+ final timestep = the final timestep value until which the per-atom quantity is to be added
<I>property/atom</I> value = varname var_value<B>0</B> var_value<B>1</B> ....
varname = name of the variable to be set
var_value<B>0</B>, var_value<B>1</B>... = values of the property to be set
@@ -84,6 +88,7 @@ set type 3 charge 0.5
set type 1*3 charge 0.5
set atom 100*200 x 0.5 y 1.0
set atom 1492 type 3
+set region reg add yes until 1000 property/atom liqOnParticle 5e-5
set property/atom Temp 273.15
</PRE>
<P><B>LIGGGHTS vs. LAMMPS Info:</B>
@@ -315,6 +320,11 @@ and the appropriate number if the property is a vector. See
<A HREF = "fix_property.html">fix property/atom</A> for details on how to define
such a per-particle property.
</P>
+<P>Keyword <I>add</I> and <I>until</I> are used to add a per-atom quantity (or property)
+in addition to the keyword <I>property/atom</I>. The <I>add</I> keyword is used to turn on
+the addition of a quantity in a region or a group
+until the timestep defined by the <I>until</I> keyword.
+</P>
<P><B>Restrictions:</B>
</P>
<P>You cannot set an atom attribute (e.g. <I>mol</I> or <I>q</I> or <I>volume</I>) if
diff --git a/doc/set.txt b/doc/set.txt
index 5ab92b7535c9eb83bce4175d7585c64f606faed1..0bcc0ff86609fa9fa3c003ff9220c65a21992889 100644
--- a/doc/set.txt
+++ b/doc/set.txt
@@ -64,6 +64,10 @@ keyword = {type} or {type/fraction} or {mol} or {x} or {y} or {z} or \
{meso_e} value = energy of SPH particles (need units)
{meso_cv} value = heat capacity of SPH particles (need units)
{meso_rho} value = density of SPH particles (need units)
+ {add} value = yes no
+ yes = add per-atom quantities to a region or a group
+ {until} value = final timestep
+ final timestep = the final timestep value until which the per-atom quantity is to be added
{property/atom} value = varname var_value[0] var_value[1] ....
varname = name of the variable to be set
var_value[0], var_value[1]... = values of the property to be set :pre
@@ -79,6 +83,7 @@ set type 3 charge 0.5
set type 1*3 charge 0.5
set atom 100*200 x 0.5 y 1.0
set atom 1492 type 3
+set region reg add yes until 1000 property/atom liqOnParticle 5e-5
set property/atom Temp 273.15 :pre
[LIGGGHTS vs. LAMMPS Info:]
@@ -310,6 +315,11 @@ and the appropriate number if the property is a vector. See
"fix property/atom"_fix_property.html for details on how to define
such a per-particle property.
+Keyword {add} and {until} are used to add a per-atom quantity (or property)
+in addition to the keyword {property/atom}. The {add} keyword is used to turn on
+the addition of a quantity in a region or a group
+until the timestep defined by the {until} keyword.
+
[Restrictions:]
You cannot set an atom attribute (e.g. {mol} or {q} or {volume}) if
diff --git a/src/atom_vec_sphere.cpp b/src/atom_vec_sphere.cpp
index a12cfae32bb42cd389f5c388d9acca3efedbadf8..811fbb22726b52e3c06d7196e4eeec126b02971a 100644
--- a/src/atom_vec_sphere.cpp
+++ b/src/atom_vec_sphere.cpp
@@ -35,6 +35,7 @@
#include "math_const.h"
#include "memory.h"
#include "error.h"
+#include "domain_wedge.h"
using namespace LAMMPS_NS;
using namespace MathConst;
@@ -238,7 +239,11 @@ int AtomVecSphere::pack_comm(int n, int *list, double *buf,
int AtomVecSphere::pack_comm_vel(int n, int *list, double *buf,
int pbc_flag, int *pbc)
{
+ if(dynamic_cast<DomainWedge*>(domain))
+ return pack_comm_vel_wedge(n,list,buf,pbc_flag,pbc);
+
int i,j,m;
+
double dx,dy,dz,dvx,dvy,dvz;
if (radvary == 0) {
@@ -280,6 +285,7 @@ int AtomVecSphere::pack_comm_vel(int n, int *list, double *buf,
buf[m++] = omega[j][2];
}
} else {
+
dvx = pbc[0]*h_rate[0] + pbc[5]*h_rate[5] + pbc[4]*h_rate[4];
dvy = pbc[1]*h_rate[1] + pbc[3]*h_rate[3];
dvz = pbc[2]*h_rate[2];
@@ -349,6 +355,7 @@ int AtomVecSphere::pack_comm_vel(int n, int *list, double *buf,
buf[m++] = omega[j][2];
}
} else {
+
dvx = pbc[0]*h_rate[0] + pbc[5]*h_rate[5] + pbc[4]*h_rate[4];
dvy = pbc[1]*h_rate[1] + pbc[3]*h_rate[3];
dvz = pbc[2]*h_rate[2];
@@ -445,6 +452,7 @@ void AtomVecSphere::unpack_comm_vel(int n, int first, double *buf)
omega[i][0] = buf[m++];
omega[i][1] = buf[m++];
omega[i][2] = buf[m++];
+
}
} else {
m = 0;
@@ -528,6 +536,7 @@ void AtomVecSphere::unpack_reverse(int n, int *list, double *buf)
m = 0;
for (i = 0; i < n; i++) {
j = list[i];
+
f[j][0] += buf[m++];
f[j][1] += buf[m++];
f[j][2] += buf[m++];
@@ -606,6 +615,9 @@ int AtomVecSphere::pack_border(int n, int *list, double *buf,
int AtomVecSphere::pack_border_vel(int n, int *list, double *buf,
int pbc_flag, int *pbc)
{
+ if(dynamic_cast<DomainWedge*>(domain))
+ return pack_border_vel_wedge(n,list,buf,pbc_flag,pbc);
+
int i,j,m;
double dx,dy,dz,dvx,dvy,dvz;
@@ -664,6 +676,7 @@ int AtomVecSphere::pack_border_vel(int n, int *list, double *buf,
dvz = pbc[2]*h_rate[2];
for (i = 0; i < n; i++) {
j = list[i];
+
buf[m++] = x[j][0] + dx;
buf[m++] = x[j][1] + dy;
buf[m++] = x[j][2] + dz;
@@ -754,6 +767,7 @@ void AtomVecSphere::unpack_border_vel(int n, int first, double *buf)
omega[i][0] = buf[m++];
omega[i][1] = buf[m++];
omega[i][2] = buf[m++];
+
}
}
diff --git a/src/atom_vec_sphere.h b/src/atom_vec_sphere.h
index a9088353bac3609c108ff6250076991cf717ffd0..a73085c2a557dab7803a9146ace0ccd3be258aa5 100644
--- a/src/atom_vec_sphere.h
+++ b/src/atom_vec_sphere.h
@@ -46,6 +46,7 @@ class AtomVecSphere : public AtomVec {
void copy(int, int, int);
int pack_comm(int, int *, double *, int, int *);
int pack_comm_vel(int, int *, double *, int, int *);
+ int pack_comm_vel_wedge(int, int *, double *, int, int *);
int pack_comm_hybrid(int, int *, double *);
void unpack_comm(int, int, double *);
void unpack_comm_vel(int, int, double *);
@@ -56,6 +57,7 @@ class AtomVecSphere : public AtomVec {
int unpack_reverse_hybrid(int, int *, double *);
int pack_border(int, int *, double *, int, int *);
int pack_border_vel(int, int *, double *, int, int *);
+ int pack_border_vel_wedge(int, int *, double *, int, int *);
int pack_border_hybrid(int, int *, double *);
void unpack_border(int, int, double *);
void unpack_border_vel(int, int, double *);
diff --git a/src/atom_vec_sphere_w.cpp b/src/atom_vec_sphere_w.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..adffed4ea73f92e46ef1d72cba8061f37829b256
--- /dev/null
+++ b/src/atom_vec_sphere_w.cpp
@@ -0,0 +1,46 @@
+/* ----------------------------------------------------------------------
+ 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 "atom_vec_sphere.h"
+#include "domain_wedge.h"
+
+#ifndef DOMAIN_WEDGE_REAL_H
+#define DOMAIN_WEDGE_REAL_H
+
+using namespace LAMMPS_NS;
+
+/* ---------------------------------------------------------------------- */
+
+int AtomVecSphere::pack_border_vel_wedge(int n, int *list, double *buf,
+ int pbc_flag, int *pbc)
+{
+ return 0;
+}
+
+/* ---------------------------------------------------------------------- */
+
+int AtomVecSphere::pack_comm_vel_wedge(int n, int *list, double *buf,
+ int pbc_flag, int *pbc)
+{
+ return 0;
+}
+
+#endif
diff --git a/src/comm.cpp b/src/comm.cpp
index 619c950b17e8767c2f0015330d76f4f32908bae5..94cae842f8008eadf88645f04d406040d2468c04 100644
--- a/src/comm.cpp
+++ b/src/comm.cpp
@@ -38,6 +38,8 @@
#include "force.h"
#include "pair.h"
#include "domain.h"
+#include "domain_wedge.h"
+#include "math_extra_liggghts.h"
#include "neighbor.h"
#include "group.h"
#include "modify.h"
@@ -63,6 +65,10 @@ enum{MULTIPLE}; // same as in ProcMap
enum{ONELEVEL,TWOLEVEL,NUMA,CUSTOM};
enum{CART,CARTREORDER,XYZ};
+// *************************************
+#include "comm_I.h"
+// *************************************
+
/* ----------------------------------------------------------------------
setup MPI and allocate buffer space
------------------------------------------------------------------------- */
@@ -546,6 +552,22 @@ void Comm::setup()
nswap = 2 * (maxneed[0]+maxneed[1]+maxneed[2]);
if (nswap > maxswap) grow_swap(nswap);
+ dw_ = dynamic_cast<DomainWedge*>(domain);
+ if(dw_)
+ {
+ ia = dw_->index_axis();
+ iphi = dw_->index_phi();
+ dw_->n1(nleft);
+ dw_->n2(nright);
+ dw_->center(c);
+ double cutghmax = MathExtraLiggghts::max(cutghost[0],cutghost[1],cutghost[2]);
+ pleft[0] = c[0] - nleft[0] * (atom->radius ? (cutghmax / 2. + neighbor->skin/2.) : cutghmax);
+ pleft[1] = c[1] - nleft[1] * (atom->radius ? (cutghmax / 2. + neighbor->skin/2.) : cutghmax);
+ pright[0] = c[0] - nright[0] * (atom->radius ? (cutghmax / 2. + neighbor->skin/2.) : cutghmax);
+ pright[1] = c[1] - nright[1] * (atom->radius ? (cutghmax / 2. + neighbor->skin/2.) : cutghmax);
+
+ }
+
// setup parameters for each exchange:
// sendproc = proc to send to at each swap
// recvproc = proc to recv from at each swap
@@ -568,6 +590,14 @@ void Comm::setup()
int iswap = 0;
for (dim = 0; dim < 3; dim++) {
for (ineed = 0; ineed < 2*maxneed[dim]; ineed++) {
+
+ if(dw_ && dim != ia && dim != iphi)
+ {
+
+ nswap--;
+ continue;
+ }
+
pbc_flag[iswap] = 0;
pbc[iswap][0] = pbc[iswap][1] = pbc[iswap][2] =
pbc[iswap][3] = pbc[iswap][4] = pbc[iswap][5] = 0;
@@ -578,7 +608,7 @@ void Comm::setup()
if (style == SINGLE) {
if (ineed < 2) slablo[iswap] = -BIG;
else slablo[iswap] = 0.5 * (sublo[dim] + subhi[dim]);
- slabhi[iswap] = sublo[dim] + cutghost[dim];
+ slabhi[iswap] = sublo[dim] + (atom->radius ? (cutghost[dim] / 2. + neighbor->skin/2.) : cutghost[dim]);
} else {
for (i = 1; i <= ntypes; i++) {
if (ineed < 2) multilo[iswap][i] = -BIG;
@@ -599,7 +629,7 @@ void Comm::setup()
sendproc[iswap] = procneigh[dim][1];
recvproc[iswap] = procneigh[dim][0];
if (style == SINGLE) {
- slablo[iswap] = subhi[dim] - cutghost[dim];
+ slablo[iswap] = subhi[dim] - (atom->radius ? (cutghost[dim] / 2. + neighbor->skin/2.) : cutghost[dim]);
if (ineed < 2) slabhi[iswap] = BIG;
else slabhi[iswap] = 0.5 * (sublo[dim] + subhi[dim]);
} else {
@@ -941,10 +971,17 @@ void Comm::borders()
smax = rmax = 0;
for (dim = 0; dim < 3; dim++) {
+
nlast = 0;
twoneed = 2*maxneed[dim];
for (ineed = 0; ineed < twoneed; ineed++) {
+ if(dw_ && dim != ia && dim != iphi)
+ {
+
+ continue;
+ }
+
// find atoms within slab boundaries lo/hi using <= and >=
// check atoms between nfirst and nlast
// for first swaps in a dim, check owned and ghost
@@ -982,11 +1019,23 @@ void Comm::borders()
if (sendflag) {
if (!bordergroup || ineed >= 2) {
if (style == SINGLE) {
- for (i = nfirst; i < nlast; i++)
- if (x[i][dim] >= lo && x[i][dim] <= hi) {
- if (nsend == maxsendlist[iswap]) grow_list(iswap,nsend);
- sendlist[iswap][nsend++] = i;
- }
+
+ if(dw_ && dim == iphi)
+ {
+ for (i = nfirst; i < nlast; i++)
+ if (decide_wedge(i,dim,lo,hi,ineed)) {
+ if (nsend == maxsendlist[iswap]) grow_list(iswap,nsend);
+ sendlist[iswap][nsend++] = i;
+ }
+ }
+ else
+ {
+ for (i = nfirst; i < nlast; i++)
+ if (decide(i,dim,lo,hi,ineed)) {
+ if (nsend == maxsendlist[iswap]) grow_list(iswap,nsend);
+ sendlist[iswap][nsend++] = i;
+ }
+ }
} else {
for (i = nfirst; i < nlast; i++) {
itype = type[i];
@@ -999,17 +1048,35 @@ void Comm::borders()
} else {
if (style == SINGLE) {
- ngroup = atom->nfirst;
- for (i = 0; i < ngroup; i++)
- if (x[i][dim] >= lo && x[i][dim] <= hi) {
- if (nsend == maxsendlist[iswap]) grow_list(iswap,nsend);
- sendlist[iswap][nsend++] = i;
- }
- for (i = atom->nlocal; i < nlast; i++)
- if (x[i][dim] >= lo && x[i][dim] <= hi) {
- if (nsend == maxsendlist[iswap]) grow_list(iswap,nsend);
- sendlist[iswap][nsend++] = i;
- }
+
+ if(dw_ && dim == iphi)
+ {
+ ngroup = atom->nfirst;
+ for (i = 0; i < ngroup; i++)
+ if (decide_wedge(i,dim,lo,hi,ineed)) {
+ if (nsend == maxsendlist[iswap]) grow_list(iswap,nsend);
+ sendlist[iswap][nsend++] = i;
+ }
+ for (i = atom->nlocal; i < nlast; i++)
+ if (decide_wedge(i,dim,lo,hi,ineed)) {
+ if (nsend == maxsendlist[iswap]) grow_list(iswap,nsend);
+ sendlist[iswap][nsend++] = i;
+ }
+ }
+ else
+ {
+ ngroup = atom->nfirst;
+ for (i = 0; i < ngroup; i++)
+ if (decide(i,dim,lo,hi,ineed)) {
+ if (nsend == maxsendlist[iswap]) grow_list(iswap,nsend);
+ sendlist[iswap][nsend++] = i;
+ }
+ for (i = atom->nlocal; i < nlast; i++)
+ if (decide(i,dim,lo,hi,ineed)) {
+ if (nsend == maxsendlist[iswap]) grow_list(iswap,nsend);
+ sendlist[iswap][nsend++] = i;
+ }
+ }
} else {
ngroup = atom->nfirst;
for (i = 0; i < ngroup; i++) {
diff --git a/src/comm.h b/src/comm.h
index 586cc5d999609a29cfecabef5b1b3377016da044..88759e5273f39f522ac833cc03ff2b4bc51595c0 100644
--- a/src/comm.h
+++ b/src/comm.h
@@ -130,6 +130,14 @@ class Comm : protected Pointers {
virtual void allocate_multi(int); // allocate multi arrays
virtual void free_swap(); // free swap arrays
virtual void free_multi(); // free multi arrays
+
+ bool decide(int i,int dim,double lo,double hi,int ineed);
+ bool decide_wedge(int i,int dim,double lo,double hi,int ineed);
+
+ class DomainWedge *dw_;
+ int ia,iphi;
+
+ double nleft[2],nright[2],pleft[2],pright[2],c[2];
};
}
diff --git a/src/comm_I.h b/src/comm_I.h
new file mode 100644
index 0000000000000000000000000000000000000000..45afbd83440c7c1755328829b7741f3c3b0d1302
--- /dev/null
+++ b/src/comm_I.h
@@ -0,0 +1,81 @@
+/* ----------------------------------------------------------------------
+ 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.
+------------------------------------------------------------------------- */
+
+#ifndef LMP_COMM_I_H
+#define LMP_COMM_I_H
+
+#include "atom.h"
+#include "domain_wedge.h"
+
+using namespace LAMMPS_NS;
+
+/* ----------------------------------------------------------------------
+ decide if border element, optimization for granular
+------------------------------------------------------------------------- */
+
+inline bool Comm::decide(int i,int dim,double lo,double hi,int ineed)
+{
+ double **x = atom->x;
+ double *radius = atom->radius;
+
+ if( ((ineed % 2 == 0) && x[i][dim] >= lo && x[i][dim] <= (hi + (radius? (radius[i]) : 0.)) ) ||
+ ((ineed % 2 == 1) && x[i][dim] >= (lo - (radius? radius[i] : 0.)) && x[i][dim] <= hi ) )
+ return true;
+
+ return false;
+}
+
+/* ----------------------------------------------------------------------
+ decide if border element for wedge case, optimization for granular
+------------------------------------------------------------------------- */
+
+inline bool Comm::decide_wedge(int i,int dim,double lo,double hi,int ineed)
+{
+ double **x = atom->x;
+ double *radius = atom->radius;
+ double coo[2],d[2];
+ coo[0] = x[i][iphi];
+ coo[1] = x[i][(iphi+1)%3];
+
+ if (ineed % 2 == 0)
+ {
+ vectorSubtract2D(coo,pleft,d);
+ if(vectorDot2D(d,nleft) >= -(radius? radius[i] : 0.))
+ {
+
+ return true;
+ }
+ }
+
+ else if (ineed % 2 == 1)
+ {
+ vectorSubtract2D(coo,pright,d);
+ if(vectorDot2D(d,nright) >= -(radius? radius[i] : 0.))
+ {
+
+ return true;
+ }
+ }
+
+ return false;
+}
+
+#endif
diff --git a/src/compute_erotate_sphere.cpp b/src/compute_erotate_sphere.cpp
index f0ba88cda6257e1bc2fdb69141b3189eda37f7c5..20e9a093c9989e681b24469cebd5bf2e340bed54 100644
--- a/src/compute_erotate_sphere.cpp
+++ b/src/compute_erotate_sphere.cpp
@@ -18,6 +18,8 @@
#include "update.h"
#include "force.h"
#include "domain.h"
+#include "modify.h"
+#include "fix_multisphere.h"
#include "group.h"
#include "error.h"
@@ -39,6 +41,8 @@ ComputeERotateSphere::ComputeERotateSphere(LAMMPS *lmp, int narg, char **arg) :
if (!atom->sphere_flag)
error->all(FLERR,"Compute erotate/sphere requires atom style sphere");
+
+ fix_ms = 0;
}
/* ---------------------------------------------------------------------- */
@@ -46,6 +50,7 @@ ComputeERotateSphere::ComputeERotateSphere(LAMMPS *lmp, int narg, char **arg) :
void ComputeERotateSphere::init()
{
pfactor = 0.5 * force->mvv2e * INERTIA;
+ fix_ms = static_cast<FixMultisphere*>(modify->find_fix_style("multisphere",0));
}
/* ---------------------------------------------------------------------- */
@@ -65,7 +70,7 @@ double ComputeERotateSphere::compute_scalar()
double erotate = 0.0;
for (int i = 0; i < nlocal; i++)
- if (mask[i] & groupbit)
+ if (mask[i] & groupbit && (!fix_ms || fix_ms->belongs_to(i) < 0))
erotate += (omega[i][0]*omega[i][0] + omega[i][1]*omega[i][1] +
omega[i][2]*omega[i][2]) * radius[i]*radius[i]*rmass[i];
diff --git a/src/compute_erotate_sphere.h b/src/compute_erotate_sphere.h
index 815a6655835a4d98a0588ab9b0741d25ed203863..3e41f304dc9afdfbf6ea57cdca2428de86c96f06 100644
--- a/src/compute_erotate_sphere.h
+++ b/src/compute_erotate_sphere.h
@@ -33,6 +33,7 @@ class ComputeERotateSphere : public Compute {
private:
double pfactor;
+ class FixMultisphere *fix_ms;
};
}
diff --git a/src/compute_ke.cpp b/src/compute_ke.cpp
index 3efd6983250a3ca922921928e77f5fe8c89ac082..65927309ff7644e49b4ddc91ddd0ea73a32c54fc 100644
--- a/src/compute_ke.cpp
+++ b/src/compute_ke.cpp
@@ -19,6 +19,8 @@
#include "domain.h"
#include "group.h"
#include "error.h"
+#include "modify.h"
+#include "fix_multisphere.h"
using namespace LAMMPS_NS;
@@ -31,6 +33,7 @@ ComputeKE::ComputeKE(LAMMPS *lmp, int narg, char **arg) :
scalar_flag = 1;
extscalar = 1;
+ fix_ms = NULL;
}
/* ---------------------------------------------------------------------- */
@@ -38,6 +41,7 @@ ComputeKE::ComputeKE(LAMMPS *lmp, int narg, char **arg) :
void ComputeKE::init()
{
pfactor = 0.5 * force->mvv2e;
+ fix_ms = static_cast<FixMultisphere*>(modify->find_fix_style("multisphere",0));
}
/* ---------------------------------------------------------------------- */
@@ -57,7 +61,7 @@ double ComputeKE::compute_scalar()
if (rmass) {
for (int i = 0; i < nlocal; i++)
- if (mask[i] & groupbit)
+ if (mask[i] & groupbit && (!fix_ms || fix_ms->belongs_to(i) < 0))
ke += rmass[i] * (v[i][0]*v[i][0] + v[i][1]*v[i][1] + v[i][2]*v[i][2]);
} else {
for (int i = 0; i < nlocal; i++)
diff --git a/src/compute_ke.h b/src/compute_ke.h
index 5dbfeb3081be3ab4de70721c150208a3abb9967e..88406c03b70ac283da1d20bc0b59a1c7382d0d28 100644
--- a/src/compute_ke.h
+++ b/src/compute_ke.h
@@ -32,6 +32,7 @@ class ComputeKE : public Compute {
private:
double pfactor;
+ class FixMultisphere *fix_ms;
};
}
diff --git a/src/compute_ke_atom.cpp b/src/compute_ke_atom.cpp
index 5709c428374be4f9e0468207a804878ddfe767ed..869d56e957048e6a2d80785665023e338765c980 100644
--- a/src/compute_ke_atom.cpp
+++ b/src/compute_ke_atom.cpp
@@ -17,6 +17,7 @@
#include "update.h"
#include "modify.h"
#include "comm.h"
+#include "fix_multisphere.h"
#include "force.h"
#include "memory.h"
#include "error.h"
@@ -35,6 +36,7 @@ ComputeKEAtom::ComputeKEAtom(LAMMPS *lmp, int narg, char **arg) :
nmax = 0;
ke = NULL;
+ fix_ms = NULL;
}
/* ---------------------------------------------------------------------- */
@@ -53,6 +55,7 @@ void ComputeKEAtom::init()
if (strcmp(modify->compute[i]->style,"ke/atom") == 0) count++;
if (count > 1 && comm->me == 0)
error->warning(FLERR,"More than one compute ke/atom");
+ fix_ms = static_cast<FixMultisphere*>(modify->find_fix_style("multisphere",0));
}
/* ---------------------------------------------------------------------- */
@@ -82,7 +85,7 @@ void ComputeKEAtom::compute_peratom()
if (rmass)
for (int i = 0; i < nlocal; i++) {
- if (mask[i] & groupbit) {
+ if (mask[i] & groupbit && (!fix_ms || fix_ms->belongs_to(i) < 0)) {
ke[i] = 0.5 * mvv2e * rmass[i] *
(v[i][0]*v[i][0] + v[i][1]*v[i][1] + v[i][2]*v[i][2]);
} else ke[i] = 0.0;
diff --git a/src/compute_ke_atom.h b/src/compute_ke_atom.h
index 3c815ed1c32eb693c4ef83223a5d5a894cca57ba..ccd91c95912c888d9d14a6d697f75b5275807a98 100644
--- a/src/compute_ke_atom.h
+++ b/src/compute_ke_atom.h
@@ -35,6 +35,7 @@ class ComputeKEAtom : public Compute {
private:
int nmax;
double *ke;
+ class FixMultisphere *fix_ms;
};
}
diff --git a/src/domain.cpp b/src/domain.cpp
index 286494f39f55eeadef66779d1e253a5a2aca09ef..d73e79f94a2d646875cd946524e2636f039c51b2 100644
--- a/src/domain.cpp
+++ b/src/domain.cpp
@@ -94,6 +94,7 @@ Domain::Domain(LAMMPS *lmp) : Pointers(lmp)
nregion = maxregion = 0;
regions = NULL;
+ is_wedge = false;
}
/* ---------------------------------------------------------------------- */
@@ -1416,41 +1417,3 @@ void Domain::box_corners()
lamda2x(corners[7],corners[7]);
}
-/* ----------------------------------------------------------------------
- domain check - not used very often, so not inlined
-------------------------------------------------------------------------- */
-
-int Domain::is_periodic_ghost(int i)
-{
- int idim;
- int nlocal = atom->nlocal;
- double *x = atom->x[i];
- double halfskin = 0.5*neighbor->skin;
-
- if(i < nlocal) return 0;
- else
- {
- for(idim = 0; idim < 3; idim++)
- if ((x[idim] < (boxlo[idim]+halfskin) || x[idim] > (boxhi[idim]-halfskin)) && periodicity[idim])
-
- return 1;
- }
- return 0;
-}
-
-/* ----------------------------------------------------------------------
- check if atom is unique on this subdomain
- used when tallying stats across owned and ghost particles
-------------------------------------------------------------------------- */
-
-bool Domain::is_owned_or_first_ghost(int i)
-{
- if(!atom->tag_enable)
- error->one(FLERR,"The current simulation setup requires atoms to have tags");
- if(0 == atom->map_style)
- error->one(FLERR,"The current simulation setup requires an 'atom_modify map' command to allocate an atom map");
-
- if(i == atom->map(atom->tag[i]))
- return true;
- return false;
-}
diff --git a/src/domain.h b/src/domain.h
index e61f6d770e22e0c15cc89c5d1839113ec8b38396..fb4f37d33312bb563c572bacd2a55c1d71b443d7 100644
--- a/src/domain.h
+++ b/src/domain.h
@@ -28,6 +28,8 @@
#include "error.h"
#include "comm.h"
#include "vector_liggghts.h"
+#include "neighbor.h"
+#include "atom.h"
#define SMALL_DMBRDR 1.0e-8
@@ -123,7 +125,7 @@ class Domain : protected Pointers {
void add_region(int, char **);
void delete_region(int, char **);
int find_region(char *);
- void set_boundary(int, char **, int);
+ virtual void set_boundary(int, char **, int);
virtual void print_box(const char *);
virtual void lamda2x(int);
@@ -134,109 +136,26 @@ class Domain : protected Pointers {
void bbox(double *, double *, double *, double *);
void box_corners();
- inline int is_in_domain(double* pos);
- inline int is_in_subdomain(double* pos);
- inline int is_in_extended_subdomain(double* pos);
- inline double dist_subbox_borders(double* pos);
+ int is_in_domain(double* pos);
+ int is_in_subdomain(double* pos);
+ int is_in_extended_subdomain(double* pos);
+ double dist_subbox_borders(double* pos);
int is_periodic_ghost(int i);
bool is_owned_or_first_ghost(int i);
+ virtual int is_in_domain_wedge(double* pos) {return 0;}
+ virtual int is_in_subdomain_wedge(double* pos) {return 0;}
+ virtual int is_in_extended_subdomain_wedge(double* pos) {return 0;}
+ virtual double dist_subbox_borders_wedge(double* pos) {return 0.;}
+ virtual int is_periodic_ghost_wedge(int i) {return 0;}
+
+ bool is_wedge;
+
private:
double small[3]; // fractions of box lengths
};
-/* ----------------------------------------------------------------------
- check if coordinate in domain, subdomain or extended subdomain
- need to test with <= and >= for domain, and < and >= for subdomain
- inlined for performance
-------------------------------------------------------------------------- */
-
-inline int Domain::is_in_domain(double* pos)
-{
- if
- (
- pos[0] >= boxlo[0] && pos[0] <= boxhi[0] &&
- pos[1] >= boxlo[1] && pos[1] <= boxhi[1] &&
- pos[2] >= boxlo[2] && pos[2] <= boxhi[2]
- ) return 1;
- return 0;
-}
-
-inline int Domain::is_in_subdomain(double* pos)
-{
- double checkhi[3];
-
- // need to test with and < and >= for subdomain
- // but need to ensure elements right at boxhi
- // are tested correctly
- if(subhi[0] == boxhi[0])
- checkhi[0] = boxhi[0] + SMALL_DMBRDR;
- else
- checkhi[0] = subhi[0];
-
- if(subhi[1] == boxhi[1])
- checkhi[1] = boxhi[1] + SMALL_DMBRDR;
- else
- checkhi[1] = subhi[1];
-
- if(subhi[2] == boxhi[2])
- checkhi[2] = boxhi[2] + SMALL_DMBRDR;
- else
- checkhi[2] = subhi[2];
-
- if ( pos[0] >= sublo[0] && pos[0] < checkhi[0] &&
- pos[1] >= sublo[1] && pos[1] < checkhi[1] &&
- pos[2] >= sublo[2] && pos[2] < checkhi[2])
- return 1;
- return 0;
-}
-
-inline int Domain::is_in_extended_subdomain(double* pos)
-{
- // called on insertion
- // yields true if particle would be in subdomain after box extension
-
- if (is_in_subdomain(pos))
- return 1;
- else if (dimension == 2)
- error->all(FLERR,"Domain::is_in_extended_subdomain() not implemented for 2d");
- else
- {
- bool flag = true;
- for(int idim = 0; idim < 3; idim++)
- {
-
- if (comm->procgrid[idim] == 1) {}
- else if(comm->myloc[idim] == comm->procgrid[idim]-1)
- flag = flag && (pos[idim] >= sublo[idim]);
- else if(comm->myloc[idim] == 0)
- flag = flag && (pos[idim] <= subhi[idim]);
-
- else
- flag = flag && (pos[idim] >= sublo[idim] && pos[idim] < subhi[idim]);
- }
- if(flag) return 1;
- return 0;
- }
- return 0;
-}
-
-inline double Domain::dist_subbox_borders(double* pos)
-{
- double deltalo[3], deltahi[3], dlo, dhi;
-
- vectorSubtract3D(sublo,pos,deltalo);
- vectorSubtract3D(subhi,pos,deltahi);
- vectorAbs3D(deltalo);
- vectorAbs3D(deltahi);
-
- dlo = vectorMin3D(deltalo);
- dhi = vectorMin3D(deltahi);
-
- if(dlo < dhi)
- return dlo;
- return dhi;
-}
+#include "domain_I.h"
}
diff --git a/src/domain_I.h b/src/domain_I.h
new file mode 100644
index 0000000000000000000000000000000000000000..54bc314aa3dacb83c55891dff7e510b6df5b27f3
--- /dev/null
+++ b/src/domain_I.h
@@ -0,0 +1,178 @@
+/* ----------------------------------------------------------------------
+ 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.
+------------------------------------------------------------------------- */
+
+#ifndef LMP_DOMAIN_I_H
+#define LMP_DOMAIN_I_H
+
+using namespace LAMMPS_NS;
+
+/* ----------------------------------------------------------------------
+ check if coordinate in domain, subdomain or extended subdomain
+ need to test with <= and >= for domain, and < and >= for subdomain
+ inlined for performance
+------------------------------------------------------------------------- */
+
+inline int Domain::is_in_domain(double* pos)
+{
+ if(is_wedge)
+ return is_in_domain_wedge(pos);
+
+ if
+ (
+ pos[0] >= boxlo[0] && pos[0] <= boxhi[0] &&
+ pos[1] >= boxlo[1] && pos[1] <= boxhi[1] &&
+ pos[2] >= boxlo[2] && pos[2] <= boxhi[2]
+ ) return 1;
+ return 0;
+}
+
+inline int Domain::is_in_subdomain(double* pos)
+{
+ if(is_wedge)
+ return is_in_subdomain_wedge(pos);
+
+ double checkhi[3];
+
+ // need to test with and < and >= for subdomain
+ // but need to ensure elements right at boxhi
+ // are tested correctly
+ if(subhi[0] == boxhi[0])
+ checkhi[0] = boxhi[0] + SMALL_DMBRDR;
+ else
+ checkhi[0] = subhi[0];
+
+ if(subhi[1] == boxhi[1])
+ checkhi[1] = boxhi[1] + SMALL_DMBRDR;
+ else
+ checkhi[1] = subhi[1];
+
+ if(subhi[2] == boxhi[2])
+ checkhi[2] = boxhi[2] + SMALL_DMBRDR;
+ else
+ checkhi[2] = subhi[2];
+
+ if ( pos[0] >= sublo[0] && pos[0] < checkhi[0] &&
+ pos[1] >= sublo[1] && pos[1] < checkhi[1] &&
+ pos[2] >= sublo[2] && pos[2] < checkhi[2])
+ return 1;
+ return 0;
+}
+
+inline int Domain::is_in_extended_subdomain(double* pos)
+{
+ if(is_wedge)
+ return is_in_extended_subdomain_wedge(pos);
+
+ // called on insertion
+ // yields true if particle would be in subdomain after box extension
+
+ if (is_in_subdomain(pos))
+ return 1;
+ else if (dimension == 2)
+ error->all(FLERR,"Domain::is_in_extended_subdomain() not implemented for 2d");
+ else
+ {
+ bool flag = true;
+ for(int idim = 0; idim < 3; idim++)
+ {
+
+ if (comm->procgrid[idim] == 1) {}
+ else if(comm->myloc[idim] == comm->procgrid[idim]-1)
+ flag = flag && (pos[idim] >= sublo[idim]);
+ else if(comm->myloc[idim] == 0)
+ flag = flag && (pos[idim] <= subhi[idim]);
+
+ else
+ flag = flag && (pos[idim] >= sublo[idim] && pos[idim] < subhi[idim]);
+ }
+ if(flag) return 1;
+ return 0;
+ }
+ return 0;
+}
+
+/* ----------------------------------------------------------------------
+ check distance from borders of subbox
+------------------------------------------------------------------------- */
+
+inline double Domain::dist_subbox_borders(double* pos)
+{
+ if(is_wedge)
+ return dist_subbox_borders_wedge(pos);
+
+ double deltalo[3], deltahi[3], dlo, dhi;
+
+ vectorSubtract3D(sublo,pos,deltalo);
+ vectorSubtract3D(subhi,pos,deltahi);
+ vectorAbs3D(deltalo);
+ vectorAbs3D(deltahi);
+
+ dlo = vectorMin3D(deltalo);
+ dhi = vectorMin3D(deltahi);
+
+ if(dlo < dhi)
+ return dlo;
+ return dhi;
+}
+
+/* ----------------------------------------------------------------------
+ domain check - not used very often, so not inlined
+------------------------------------------------------------------------- */
+
+inline int Domain::is_periodic_ghost(int i)
+{
+ if(is_wedge)
+ return is_periodic_ghost_wedge(i);
+
+ int idim;
+ int nlocal = atom->nlocal;
+ double *x = atom->x[i];
+ double halfskin = 0.5*neighbor->skin;
+
+ if(i < nlocal) return 0;
+ else
+ {
+ for(idim = 0; idim < 3; idim++)
+ if ((x[idim] < (boxlo[idim]+halfskin) || x[idim] > (boxhi[idim]-halfskin)) && periodicity[idim])
+
+ return 1;
+ }
+ return 0;
+}
+
+/* ----------------------------------------------------------------------
+ check if atom is the true representation of the particle on this subdomain
+ used when tallying stats across owned and ghost particles
+------------------------------------------------------------------------- */
+
+inline bool Domain::is_owned_or_first_ghost(int i)
+{
+ if(!atom->tag_enable)
+ error->one(FLERR,"The current simulation setup requires atoms to have tags");
+ if(0 == atom->map_style)
+ error->one(FLERR,"The current simulation setup requires an 'atom_modify map' command to allocate an atom map");
+
+ if(i == atom->map(atom->tag[i]))
+ return true;
+ return false;
+}
+
+#endif
diff --git a/src/domain_wedge_dummy.h b/src/domain_wedge_dummy.h
index 8a6bee98aaa9c9bbcd872cb9ae86cb241ddc606c..6026d614a5c5c93294edd418cc4042246ba5ee35 100644
--- a/src/domain_wedge_dummy.h
+++ b/src/domain_wedge_dummy.h
@@ -40,6 +40,20 @@ class DomainWedge : public Domain
DomainWedge(class LAMMPS *lmp) : Domain(lmp) {};
void set_domain(class RegWedge *rw) {}
+ inline int index_axis()
+ { return 0; }
+
+ inline int index_phi()
+ { return 0; }
+
+ inline void n1(double *_n1)
+ { }
+
+ inline void n2(double *_n2)
+ { }
+
+ inline void center(double *_c)
+ { }
};
}
diff --git a/src/fix.h b/src/fix.h
index ca003993fa29f48b27197bd0c910f7b295244c0e..4d3b27a3caeec245e078b0183d1d46f8c4d857af 100644
--- a/src/fix.h
+++ b/src/fix.h
@@ -101,7 +101,9 @@ class Fix : protected Pointers {
virtual void post_create_pre_restart() {}
virtual void post_create() {}
virtual void pre_delete(bool) {}
- virtual void box_extent(double &xlo,double &xhi,double &ylo,double &yhi,double &zlo,double &zhi) {}
+ virtual void box_extent(double &xlo,double &xhi,double &ylo,double &yhi,double &zlo,double &zhi) {
+ UNUSED(xlo); UNUSED(xhi); UNUSED(ylo); UNUSED(yhi); UNUSED(zlo); UNUSED(zhi);
+ }
virtual void init() {}
virtual void init_list(int, class NeighList *) {}
virtual void setup(int) {}
diff --git a/src/fix_cfd_coupling_force.cpp b/src/fix_cfd_coupling_force.cpp
index eefec9195430b6fc0091cafc191fc536584a060d..6b5326fa2025d840c7f431cea8264354c5944758 100644
--- a/src/fix_cfd_coupling_force.cpp
+++ b/src/fix_cfd_coupling_force.cpp
@@ -78,7 +78,7 @@ FixCfdCouplingForce::FixCfdCouplingForce(LAMMPS *lmp, int narg, char **arg) : Fi
error->fix_error(FLERR,this,"expecting 'yes' or 'no' after 'transfer_type'");
iarg++;
hasargs = true;
- } else
+ } else if (strcmp(this->style,"couple/cfd/force") == 0)
error->fix_error(FLERR,this,"unknown keyword");
}
diff --git a/src/fix_cfd_coupling_force_implicit.cpp b/src/fix_cfd_coupling_force_implicit.cpp
index 49e41c3ef3a60724949711183ed909c38de64739..1384f4fb2aca96f01b130b06c339472c59182497 100644
--- a/src/fix_cfd_coupling_force_implicit.cpp
+++ b/src/fix_cfd_coupling_force_implicit.cpp
@@ -22,6 +22,7 @@
#include "string.h"
#include "stdlib.h"
#include "atom.h"
+#include "force.h"
#include "update.h"
#include "respa.h"
#include "error.h"
@@ -41,10 +42,33 @@ using namespace FixConst;
FixCfdCouplingForceImplicit::FixCfdCouplingForceImplicit(LAMMPS *lmp, int narg, char **arg) :
FixCfdCouplingForce(lmp,narg,arg),
+ useCN_(false),
+ CNalpha_(0.0),
fix_Ksl_(0),
fix_uf_(0)
{
+ int iarg = 3;
+
+ bool hasargs = true;
+ while(iarg < narg && hasargs)
+ {
+ hasargs = false;
+
+ if(strcmp(arg[iarg],"CrankNicolson") == 0) {
+ if(narg < iarg+2)
+ error->fix_error(FLERR,this,"not enough arguments for 'CrankNicholson'");
+ iarg++;
+ useCN_ = true;
+ CNalpha_ = atof(arg[iarg]);
+ fprintf(screen,"cfd_coupling_foce_implicit will use Crank-Nicholson scheme with %f\n", CNalpha_);
+ iarg++;
+ hasargs = true;
+ }
+ }
+
+ nevery = 1;
+
}
/* ---------------------------------------------------------------------- */
@@ -54,6 +78,15 @@ FixCfdCouplingForceImplicit::~FixCfdCouplingForceImplicit()
}
+/* ---------------------------------------------------------------------- */
+int FixCfdCouplingForceImplicit::setmask()
+{
+ int mask = 0;
+ mask |= POST_FORCE;
+ mask |= END_OF_STEP;
+ return mask;
+}
+
/* ---------------------------------------------------------------------- */
void FixCfdCouplingForceImplicit::post_create()
@@ -64,7 +97,7 @@ void FixCfdCouplingForceImplicit::post_create()
// register Ksl
if(!fix_Ksl_)
{
- char* fixarg[9];
+ char* fixarg[9];
fixarg[0]="Ksl";
fixarg[1]="all";
fixarg[2]="property/atom";
@@ -113,12 +146,15 @@ void FixCfdCouplingForceImplicit::init()
// values to come from OF
fix_coupling_->add_pull_property("Ksl","scalar-atom");
fix_coupling_->add_pull_property("uf","vector-atom");
+
+ deltaT_ = 0.5 * update->dt * force->ftm2v;
}
/* ---------------------------------------------------------------------- */
void FixCfdCouplingForceImplicit::post_force(int vflag)
{
+
double **x = atom->x;
double **v = atom->v;
double **f = atom->f;
@@ -126,7 +162,7 @@ void FixCfdCouplingForceImplicit::post_force(int vflag)
int nlocal = atom->nlocal;
double *Ksl = fix_Ksl_->vector_atom;
double **uf = fix_uf_->array_atom;
- double **dragforce = fix_dragforce_->array_atom;
+ double **dragforce = fix_dragforce_->array_atom;
double frc[3];
vectorZeroize3D(dragforce_total);
@@ -137,16 +173,86 @@ void FixCfdCouplingForceImplicit::post_force(int vflag)
if (mask[i] & groupbit)
{
// calc force
- vectorSubtract3D(uf[i],v[i],frc);
- vectorScalarMult3D(frc,Ksl[i]);
+ if(!useCN_) //calculate drag force and add if not using Crank-Nicolson
+ {
+ vectorSubtract3D(uf[i],v[i],frc);
+ vectorScalarMult3D(frc,Ksl[i]);
+
+ vectorAdd3D(f[i],frc,f[i]);
+ vectorAdd3D(dragforce_total,frc,dragforce_total);
+ }
+
+ // add other force
+ vectorAdd3D(f[i],dragforce[i],f[i]);
+
+ // add up forces for post-proc
+ vectorAdd3D(dragforce_total,dragforce[i],dragforce_total);
+ }
+ }
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixCfdCouplingForceImplicit::end_of_step()
+{
+
+ if(!useCN_) return; //return if CN not used
+
+ double **x = atom->x;
+ double **v = atom->v;
+ double **f = atom->f;
+ double *rmass = atom->rmass;
+ double *mass = atom->mass;
+ int *type = atom->type;
+
+ int *mask = atom->mask;
+ int nlocal = atom->nlocal;
+ double *Ksl = fix_Ksl_->vector_atom;
+ double **uf = fix_uf_->array_atom;
+ double **dragforce = fix_dragforce_->array_atom;
+ double KslMDeltaT, deltaU;
+ double vN32[3];
+ double frc[3];
- // add force
+ vectorZeroize3D(dragforce_total);
+
+ // add dragforce to force vector
+ for (int i = 0; i < nlocal; i++)
+ {
+ if (mask[i] & groupbit)
+ {
+ if (rmass) KslMDeltaT = Ksl[i]/rmass[i]*deltaT_;
+ else KslMDeltaT = Ksl[i]/mass[type[i]]*deltaT_;
+
+ for(int dirI=0;dirI<3;dirI++)
+ {
+ //calculate new velocity
+ vN32[dirI] = ( v[i][dirI]
+ + KslMDeltaT
+ *( uf[i][dirI]
+ - (1.0-CNalpha_)*v[i][dirI]
+ )
+ )
+ /
+ (1.0+KslMDeltaT*CNalpha_);
+
+ //calculate velocity difference and force
+ deltaU = uf[i][dirI]
+ - (
+ (1.0-CNalpha_)*v[i][dirI]
+ + CNalpha_ *vN32[dirI]
+ );
+ frc[dirI] = Ksl[i] * deltaU; //force required for the next time step
+
+ //update the particle velocity
+ v[i][dirI] += KslMDeltaT/2.0 * deltaU; //update velocity for a half step!
+ }
+
+ // add force
vectorAdd3D(f[i],frc,f[i]);
- vectorAdd3D(f[i],dragforce[i],f[i]);
-
- // add up forces for post-proc
+
+ // add up forces for post-proc
vectorAdd3D(dragforce_total,frc,dragforce_total);
- vectorAdd3D(dragforce_total,dragforce[i],dragforce_total);
- }
+ }
}
}
diff --git a/src/fix_cfd_coupling_force_implicit.h b/src/fix_cfd_coupling_force_implicit.h
index e4918fc89e50eb7c967b07368a09370cef822286..7beddfab77839bfccfca26107357b41ba635c163 100644
--- a/src/fix_cfd_coupling_force_implicit.h
+++ b/src/fix_cfd_coupling_force_implicit.h
@@ -39,10 +39,15 @@ class FixCfdCouplingForceImplicit : public FixCfdCouplingForce {
void post_create();
void pre_delete(bool unfixflag);
+ int setmask();
virtual void init();
void post_force(int);
+ void end_of_step();
protected:
+ double deltaT_;
+ bool useCN_;
+ double CNalpha_;
class FixPropertyAtom* fix_Ksl_;
class FixPropertyAtom* fix_uf_;
};
diff --git a/src/fix_contact_history.cpp b/src/fix_contact_history.cpp
index 2c1374032273710b004c06bf0627b0f1d4ebcfe7..301f2f9bdbc3a0bf005f9a074cfaac613297f40c 100644
--- a/src/fix_contact_history.cpp
+++ b/src/fix_contact_history.cpp
@@ -214,7 +214,7 @@ void FixContactHistory::init()
only invoke pre_exchange() if neigh list stores more current history info
than npartner/partner arrays in this fix
that will only be case if pair->compute() has been invoked since
- upate of npartner/npartner
+ update of npartner/npartner
this logic avoids 2 problems:
run 100; write_restart; run 100
setup_pre_exchange is called twice (by write_restart and 2nd run setup)
diff --git a/src/fix_heat_gran.cpp b/src/fix_heat_gran.cpp
index f43fd78d673e40ed2a34306e5ffaa1144808bd96..ceb995e92fee75a0e8b647360790aa93b448fc41 100644
--- a/src/fix_heat_gran.cpp
+++ b/src/fix_heat_gran.cpp
@@ -55,7 +55,6 @@ FixHeatGran::FixHeatGran(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg
peratom_flag = 1;
size_peratom_cols = 0;
peratom_freq = 1;
- time_depend = 1;
scalar_flag = 1;
global_freq = 1;
diff --git a/src/fix_insert.cpp b/src/fix_insert.cpp
index 332a934fa7ccd522feae88ac6e9de0eef52fc83c..9dce5dde4e8713d994c262f50dbd0612f02a4631 100644
--- a/src/fix_insert.cpp
+++ b/src/fix_insert.cpp
@@ -260,6 +260,7 @@ FixInsert::FixInsert(LAMMPS *lmp, int narg, char **arg) :
vector_flag = 1;
size_vector = 2;
+ global_freq = 1;
print_stats_start_flag = 1;
@@ -298,11 +299,11 @@ void FixInsert::setup(int vflag)
// calc last step of insertion
if(ninsert_exists)
{
- if(ninsert < ninsert_per)
+ if(ninsert <= ninsert_per)
final_ins_step = first_ins_step;
else
final_ins_step = first_ins_step +
- static_cast<int>(static_cast<double>(ninsert)/ninsert_per * static_cast<double>(insert_every));
+ static_cast<int>(static_cast<double>(ninsert)/ninsert_per) * static_cast<double>(insert_every);
if(final_ins_step < 0)
error->fix_error(FLERR,this,"Particle insertion: Overflow - need too long for particle insertion. "
diff --git a/src/fix_massflow_mesh.cpp b/src/fix_massflow_mesh.cpp
index 2071f6c6005e36b147544e43f8135ce79a873758..83c044305d4c1d1cc2c002314a66f5bb6fe02658 100644
--- a/src/fix_massflow_mesh.cpp
+++ b/src/fix_massflow_mesh.cpp
@@ -476,7 +476,7 @@ void FixMassflowMesh::pre_exchange()
void FixMassflowMesh::write_restart(FILE *fp)
{
int n = 0;
- double list[4];
+ double list[6];
list[n++] = mass_;
list[n++] = t_count_;
list[n++] = mass_last_;
diff --git a/src/fix_mesh.cpp b/src/fix_mesh.cpp
index cb64b3fc8a6f7869dc54d595c9f67c1d2a07dc9e..05f87c199cf8fecec21cd80d00c6fc737a83f341 100644
--- a/src/fix_mesh.cpp
+++ b/src/fix_mesh.cpp
@@ -115,7 +115,6 @@ FixMesh::FixMesh(LAMMPS *lmp, int narg, char **arg)
}
// construct a mesh - can be surface or volume mesh
-
// just create object and return if reading data from restart file
if(modify->have_restart_data(this)) create_mesh_restart();
@@ -304,7 +303,6 @@ void FixMesh::initialSetup()
void FixMesh::pre_exchange()
{
- // flag parallel operations on this step
pOpFlag_ = true;
}
@@ -315,7 +313,6 @@ void FixMesh::pre_exchange()
void FixMesh::pre_force(int vflag)
{
-
// case re-neigh step
if(pOpFlag_)
{
@@ -346,6 +343,7 @@ void FixMesh::pre_force(int vflag)
void FixMesh::final_integrate()
{
+
mesh_->reverseComm();
}
diff --git a/src/fix_mesh_surface_stress.cpp b/src/fix_mesh_surface_stress.cpp
index 3f10759abb6449402d7e81d4e562d9e407d326ff..a75f53d9b1113cce153ce37e44d2a254d2a1a5c6 100644
--- a/src/fix_mesh_surface_stress.cpp
+++ b/src/fix_mesh_surface_stress.cpp
@@ -37,11 +37,11 @@
#include "fix_property_global.h"
#include "fix_gravity.h"
-#define EPSILON 0.001
-
using namespace LAMMPS_NS;
using namespace FixConst;
+#define EPSILON 0.0001
+
/* ---------------------------------------------------------------------- */
FixMeshSurfaceStress::FixMeshSurfaceStress(LAMMPS *lmp, int narg, char **arg)
@@ -228,7 +228,7 @@ void FixMeshSurfaceStress::final_integrate()
void FixMeshSurfaceStress::add_particle_contribution(int ip,double *frc,
double *delta,int iTri,double *v_wall)
{
- double E,c[3],v_rel[3],cmag,vmag,cos_gamma,sin_gamma,sin_2gamma,tan_gamma;
+ double E,c[3],v_rel[3],cmag,v_rel_mag,cos_gamma,sin_gamma,sin_2gamma,tan_gamma;
double contactPoint[3],surfNorm[3], tmp[3], tmp2[3];
// do not include if not in fix group
@@ -265,18 +265,23 @@ void FixMeshSurfaceStress::add_particle_contribution(int ip,double *frc,
vectorSubtract3D(v,v_wall,v_rel);
if(vectorDot3D(c,v_rel) < 0.) return;
- vmag = vectorMag3D(v_rel);
+ v_rel_mag = vectorMag3D(v_rel);
// get surface normal
+
triMesh()->surfaceNorm(iTri,surfNorm);
- sin_gamma = MathExtraLiggghts::abs(vectorDot3D(v_rel,surfNorm)) / (vmag);
-
+ // return if no relative velocity
+ if(0.0000001 > v_rel_mag)
+ return;
+
+ sin_gamma = MathExtraLiggghts::abs(vectorDot3D(v_rel,surfNorm)) / (v_rel_mag);
+ cos_gamma = MathExtraLiggghts::abs(vectorCrossMag3D(v_rel,surfNorm)) / (v_rel_mag);
+
+ if(cos_gamma > 1.) cos_gamma = 1.;
if(sin_gamma > 1.) sin_gamma = 1.;
- if(sin_gamma < -1.) sin_gamma = -1.;
- cos_gamma = sqrt(1. - sin_gamma * sin_gamma);
- if(cos_gamma > EPSILON || (sin_gamma < 3. * cos_gamma))
+ if(cos_gamma < EPSILON || 3.*sin_gamma > cos_gamma)
{
E = 0.33333 * cos_gamma * cos_gamma;
@@ -287,7 +292,7 @@ void FixMeshSurfaceStress::add_particle_contribution(int ip,double *frc,
E = sin_2gamma - 3. * sin_gamma * sin_gamma;
}
- E *= 2.*k_finnie_[atomTypeWall()-1][atom->type[ip]-1] * vmag * vectorMag3D(frc);
+ E *= 2.*k_finnie_[atomTypeWall()-1][atom->type[ip]-1] * v_rel_mag * vectorMag3D(frc);
wear_step(iTri) += E / triMesh()->areaElem(iTri);
}
diff --git a/src/fix_mesh_surface_stress_servo.cpp b/src/fix_mesh_surface_stress_servo.cpp
index 7259969919f495920041df1d0d7ad9ab87cc05ff..648564c267502687915fe610fac5ce1c4007efa9 100644
--- a/src/fix_mesh_surface_stress_servo.cpp
+++ b/src/fix_mesh_surface_stress_servo.cpp
@@ -143,7 +143,7 @@ FixMeshSurfaceStressServo::FixMeshSurfaceStressServo(LAMMPS *lmp, int narg, char
} else {
set_point_ = -force->numeric(arg[iarg_]); // the resultant force/torque/shear acts in opposite direction --> negative value
if (set_point_ == 0.) error->fix_error(FLERR,this,"'target_val' (desired force/torque) has to be != 0.0");
- set_point_inv_ = 1./abs(set_point_);
+ set_point_inv_ = 1./fabs(set_point_);
sp_style_ = CONSTANT;
}
iarg_++;
@@ -427,7 +427,7 @@ void FixMeshSurfaceStressServo::final_integrate()
set_point_ = -input->variable->compute_equal(sp_var_);
if (set_point_ == 0.) error->fix_error(FLERR,this,"Set point (desired force/torque/shear) has to be != 0.0");
- set_point_inv_ = 1./abs(set_point_);
+ set_point_inv_ = 1./fabs(set_point_);
modify->addstep_compute(update->ntimestep + 1);
@@ -449,12 +449,12 @@ void FixMeshSurfaceStressServo::final_integrate()
// cruise mode
test_output = ctrl_output_max_;
} else {
- if (abs(err_) <= e_low) {
+ if (fabs(err_) <= e_low) {
test_output = tmp_scale*ctrl_output_min_;
- } else if(abs(err_) >= e_high) {
+ } else if(fabs(err_) >= e_high) {
test_output = ctrl_output_min_;
} else { // linear interpolation
- test_output = tmp_scale*ctrl_output_min_ + ((1-tmp_scale)*ctrl_output_min_) * (abs(err_)-e_low)/(e_high-e_low);
+ test_output = tmp_scale*ctrl_output_min_ + ((1-tmp_scale)*ctrl_output_min_) * (fabs(err_)-e_low)/(e_high-e_low);
}
}
@@ -469,7 +469,7 @@ void FixMeshSurfaceStressServo::final_integrate()
set_point_ = -input->variable->compute_equal(sp_var_);
if (set_point_ == 0.) error->fix_error(FLERR,this,"Set point (desired force/torque/shear) has to be != 0.0");
- set_point_inv_ = 1./abs(set_point_);
+ set_point_inv_ = 1./fabs(set_point_);
modify->addstep_compute(update->ntimestep + 1);
@@ -514,7 +514,7 @@ void FixMeshSurfaceStressServo::limit_vel()
{
double vmag, factor, maxOutput;
- vmag = abs(*control_output_);
+ vmag = fabs(*control_output_);
// saturation of the velocity
int totNumContacts = fix_mesh_neighlist_->getTotalNumContacts();
@@ -612,7 +612,7 @@ int FixMeshSurfaceStressServo::modify_param(int narg, char **arg)
} else {
set_point_ = -force->numeric(arg[1]); // the resultant force/torque/shear acts in opposite direction --> negative value
if (set_point_ == 0.) error->fix_error(FLERR,this,"'target_val' (desired force/torque) has to be != 0.0");
- set_point_inv_ = 1./abs(set_point_);
+ set_point_inv_ = 1./fabs(set_point_);
sp_style_ = CONSTANT;
}
diff --git a/src/fix_property_global.cpp b/src/fix_property_global.cpp
index d10617192e4eb61e16bb1855fee3e846a8e85150..4edb8c0042e99d44a9f204b7f74cbc82c6acf354 100644
--- a/src/fix_property_global.cpp
+++ b/src/fix_property_global.cpp
@@ -375,6 +375,7 @@ int FixPropertyGlobal::modify_param(int narg, char **arg)
filename = new char[strlen(arg[1])+1];
strcpy(filename,arg[1]);
+ grpname = new char[strlen(group->names[igroup])+1];
strcpy(grpname,group->names[igroup]);
return 2;
}
diff --git a/src/lammps.cpp b/src/lammps.cpp
index 7d95660809b9277e77a5b4247fbf56b5e7cd2ce0..e8500e23c0d0547731d6b55046b076eb565b4c3b 100644
--- a/src/lammps.cpp
+++ b/src/lammps.cpp
@@ -100,7 +100,7 @@ LAMMPS::LAMMPS(int narg, char **arg, MPI_Comm communicator)
universe->add_world(arg[iarg]);
iarg++;
}
- } else if (strcmp(arg[iarg],"-uid") == 0) {
+ } else if (strcmp(arg[iarg],"-uid") == 0) {
if (iarg+2 > narg)
error->universe_all(FLERR,"Invalid command-line argument");
universe->id(arg[iarg+1]);
@@ -486,6 +486,7 @@ void LAMMPS::create()
if (cuda) domain = new DomainCuda(this);
else if (wedgeflag) domain = new DomainWedge(this);
else domain = new Domain(this);
+ domain->is_wedge = wedgeflag;
group = new Group(this);
force = new Force(this); // must be after group, to create temperature
diff --git a/src/math_extra_liggghts.h b/src/math_extra_liggghts.h
index 3104dbea4b86884e6dd12ada58ab63daf5cf0cb9..9fb6adc7db99eea1ecdc88c56cbf52ce793b05f4 100644
--- a/src/math_extra_liggghts.h
+++ b/src/math_extra_liggghts.h
@@ -22,6 +22,7 @@
#ifndef LMP_MATH_EXTRA_LIGGGHTS_H
#define LMP_MATH_EXTRA_LIGGGHTS_H
+#include "pointers.h"
#include "math.h"
#include "stdio.h"
#include "string.h"
@@ -70,8 +71,8 @@ namespace MathExtraLiggghts {
inline void vec_from_quat(const double *q, double *v);
inline void vec_quat_rotate(double *vec, double *quat, double *result);
inline void vec_quat_rotate(double *vec, double *quat);
- inline void vec_quat_rotate(int *vec, double *quat) {}
- inline void vec_quat_rotate(bool *vec, double *quat) {}
+ inline void vec_quat_rotate(int *vec, double *quat) { UNUSED(vec); UNUSED(quat); }
+ inline void vec_quat_rotate(bool *vec, double *quat) { UNUSED(vec); UNUSED(quat); }
inline void quat_diff(double *q_new, double *q_old, double *q_diff);
inline void angmom_from_omega(double *w,
double *ex, double *ey, double *ez,
@@ -105,6 +106,7 @@ Matrix determinant
double MathExtraLiggghts::mdet(const double m[3][3],LAMMPS_NS::Error *error)
{
+ UNUSED(error);
return ( -m[0][2]*m[1][1]*m[2][0] + m[0][1]*m[1][2]*m[2][0] + m[0][2]*m[1][0]*m[2][1] - m[0][0]*m[1][2]*m[2][1] - m[0][1]*m[1][0]*m[2][2] + m[0][0]*m[1][1]*m[2][2] );
}
@@ -304,6 +306,7 @@ void MathExtraLiggghts::local_coosys_to_cartesian(double *global,double *local,
void MathExtraLiggghts::cartesian_coosys_to_local(double *local,double *global, double *ex_local, double *ey_local, double *ez_local,LAMMPS_NS::Error *error)
{
+ UNUSED(error);
double M[3][3] = {{0.,0.,0.},{0.,0.,0.},{0.,0.,0.}};
double Mt[3][3] = {{0.,0.,0.},{0.,0.,0.},{0.,0.,0.}};
@@ -512,6 +515,7 @@ void MathExtraLiggghts::calcBaryTriCoords(double *ap, double **edgeVec, double *
void MathExtraLiggghts::calcBaryTriCoords(double *ap, double *edgeVec0, double *edgeVec1, double *edgeVec2,
double *edgeLen, double *bary)
{
+ UNUSED(edgeVec1);
double a = LAMMPS_NS::vectorDot3D(ap,edgeVec0);
double b = LAMMPS_NS::vectorDot3D(ap,edgeVec2);
diff --git a/src/modify.cpp b/src/modify.cpp
index 31ca3cc3fa1352dd7084a3a791bb4227a77a82c0..d9d619cef84e2bf9933cd97f76ac204d13785f21 100644
--- a/src/modify.cpp
+++ b/src/modify.cpp
@@ -735,6 +735,8 @@ void Modify::add_fix(int narg, char **arg, char *suffix)
// if yes, pass state info to the Fix so it can reset itself
for (int i = 0; i < nfix_restart_global; i++)
+ {
+
if (strcmp(id_restart_global[i],fix[ifix]->id) == 0 &&
strcmp(style_restart_global[i],fix[ifix]->style) == 0) {
fix[ifix]->restart(state_restart_global[i]);
@@ -746,6 +748,7 @@ void Modify::add_fix(int narg, char **arg, char *suffix)
if (logfile) fprintf(logfile,str,fix[ifix]->id,fix[ifix]->style);
}
}
+ }
// check if Fix is in restart_peratom list
// if yes, loop over atoms so they can extract info from atom->extra array
diff --git a/src/modify.h b/src/modify.h
index c4852e7fc2c77aefd47f0a58eb964ad9e073efd2..f08099a2dfd2c52da733ab8103c9409902f6a406 100644
--- a/src/modify.h
+++ b/src/modify.h
@@ -107,6 +107,7 @@ class Modify : protected Pointers {
class Fix* find_fix_style(const char *style, int rank);
class Fix* find_fix_style_strict(const char *style, int rank);
int n_fixes_style(const char *style);
+ int n_computes_style(const char *style);
int n_fixes_style_strict(const char *style);
bool i_am_first_of_style(class Fix *fix_to_check);
int index_first_fix_of_style(const char *style);
diff --git a/src/modify_liggghts.cpp b/src/modify_liggghts.cpp
index b03a08610b9903ff96ce964f6a9fcdad8333f27e..4b014c617707af120fdd3d837d23eda89b66fe7d 100644
--- a/src/modify_liggghts.cpp
+++ b/src/modify_liggghts.cpp
@@ -124,7 +124,7 @@ Fix* Modify::find_fix_id_style(const char *id,const char* style)
}
/* ----------------------------------------------------------------------
- return number of fixes with requested style
+ return number of fixes/computes with requested style
------------------------------------------------------------------------- */
int Modify::n_fixes_style(const char *style)
@@ -141,6 +141,20 @@ int Modify::n_fixes_style(const char *style)
return n_fixes;
}
+int Modify::n_computes_style(const char *style)
+{
+ int n_computes,len;
+
+ n_computes = 0;
+ len = strlen(style);
+
+ for(int icompute = 0; icompute < ncompute; icompute++)
+ if(strncmp(compute[icompute]->style,style,len) == 0)
+ n_computes++;
+
+ return n_computes;
+}
+
int Modify::n_fixes_style_strict(const char *style)
{
int n_fixes,len;
diff --git a/src/multi_node_mesh.h b/src/multi_node_mesh.h
index 67a4040eb4af19f007485ea87871a092e0fc0b90..5d3e0e2e194fa165e86c16bfb49273109ce2b2d7 100644
--- a/src/multi_node_mesh.h
+++ b/src/multi_node_mesh.h
@@ -142,10 +142,10 @@ namespace LAMMPS_NS
bool nodesAreEqual(int iSurf, int iNode, int jSurf, int jNode);
bool nodesAreEqual(double *nodeToCheck1,double *nodeToCheck2);
- // returns true if surfaces share a node
- // called with local index
- // iNode, jNode return indices of first shared node
- bool shareNode(int iElem, int jElem, int &iNode, int &jNode);
+ // returns true if surfaces share 2 nodes
+ // called with local element indices
+ // returns indices of nodes with iNode1 < iNode2
+ bool share2Nodes(int iElem, int jElem, int &iNode1, int &jNode1, int &iNode2, int &jNode2);
// returns number of shared nodes
// called with local index
diff --git a/src/multi_node_mesh_I.h b/src/multi_node_mesh_I.h
index 3d1be6a62854b54cf8a149c6e3a2af89996745c2..d995e3c415d43c1bff22418c54746dfc6d66ebf1 100644
--- a/src/multi_node_mesh_I.h
+++ b/src/multi_node_mesh_I.h
@@ -228,19 +228,22 @@
}
/* ----------------------------------------------------------------------
- return the lowest iNode/jNode combination that is shared
+ return if elemens share node, returns lowest iNode and corresponding jNode
------------------------------------------------------------------------- */
template<int NUM_NODES>
- bool MultiNodeMesh<NUM_NODES>::shareNode(int iElem, int jElem, int &iNode, int &jNode)
+ bool MultiNodeMesh<NUM_NODES>::share2Nodes(int iElem, int jElem,
+ int &iNode1, int &jNode1, int &iNode2, int &jNode2)
{
// broad phase
double dist[3], radsum;
+ int nShared = 0;
vectorSubtract3D(center_(iElem),center_(jElem),dist);
radsum = rBound_(iElem) + rBound_(jElem);
+
if(vectorMag3DSquared(dist) > radsum*radsum)
{
- iNode = -1; jNode = -1;
+ iNode1 = jNode1 = iNode2 = jNode2 = -1;
return false;
}
@@ -249,13 +252,24 @@
for(int i=0;i<NUM_NODES;i++){
for(int j=0;j<NUM_NODES;j++){
if(MultiNodeMesh<NUM_NODES>::nodesAreEqual(iElem,i,jElem,j)){
- iNode = i; jNode = j;
-
- return true;
+ if(0 == nShared)
+ {
+ iNode1 = i;
+ jNode1 = j;
+ }
+ else
+ {
+ iNode2 = i;
+ jNode2 = j;
+
+ return true;
+ }
+ nShared++;
}
}
}
- iNode = -1; jNode = -1;
+
+ iNode1 = jNode1 = iNode2 = jNode2 = -1;
return false;
}
diff --git a/src/neigh_gran.cpp b/src/neigh_gran.cpp
index 78f69956b0610df3d0b6b41452c8fc724d22430b..8de54ae5bf876a232bb1e07c8ea5c71f7de16a17 100644
--- a/src/neigh_gran.cpp
+++ b/src/neigh_gran.cpp
@@ -387,7 +387,7 @@ void Neighbor::granular_bin_no_newton(NeighList *list)
rsq = delx*delx + dely*dely + delz*delz;
radsum = radi + radius[j];
cutsq = (radsum+skin) * (radsum+skin);
-
+
if (rsq <= cutsq) {
neighptr[n] = j;
diff --git a/src/neighbor.cpp b/src/neighbor.cpp
index ba0f1fa6c1a28cbbdf89d9da6de1ea24d8d1b4f2..8fd321c93479e04e2a77a5f72998a915916a6af9 100644
--- a/src/neighbor.cpp
+++ b/src/neighbor.cpp
@@ -1225,7 +1225,10 @@ int Neighbor::check_distance()
delta = 0.5 * (skin - (delta1+delta2));
deltasq = delta*delta;
}
- } else deltasq = triggersq;
+ } else {
+ deltasq = triggersq;
+ delta = sqrt(deltasq);
+ }
double **x = atom->x;
double *radius = atom->radius;
@@ -1255,7 +1258,7 @@ int Neighbor::check_distance()
delr = radius[i] - rhold[i];
rsq = delx*delx + dely*dely + delz*delz;
delrsq = delr*delr;
- if (delrsq > deltasq || rsq > deltasq - 2.*delr*delta + delr*delr)
+ if (delrsq > deltasq || rsq > deltasq - 2.*delr*delta + delrsq)
flag = 1;
}
diff --git a/src/pair_gran.cpp b/src/pair_gran.cpp
index 323427b8c1a9d8f7a9055765c09c96aedffd032c..6e2c7eb4d9c1bca6980ddb69bb9b09f8ea8e2867 100644
--- a/src/pair_gran.cpp
+++ b/src/pair_gran.cpp
@@ -192,7 +192,7 @@ void PairGran::init_style()
// error and warning checks
- if(0 != neighbor->delay)
+ if(0 == comm->me && 0 != neighbor->delay)
error->warning(FLERR,"It is heavily recommended to use 'neigh_modify delay 0' with granular pair styles");
if(strcmp(update->unit_style,"metal") ==0 || strcmp(update->unit_style,"real") == 0)
diff --git a/src/pair_gran_hooke.cpp b/src/pair_gran_hooke.cpp
index 4594f0039845501dddf6646d99ffbda4863cb828..249a536f11a868c7e1d22a2c7d5cb50e2f9dfa8b 100644
--- a/src/pair_gran_hooke.cpp
+++ b/src/pair_gran_hooke.cpp
@@ -29,6 +29,7 @@
#include "pair_gran_hooke.h"
#include "atom.h"
#include "force.h"
+#include "update.h"
#include "neigh_list.h"
#include "error.h"
#include "vector_liggghts.h"
@@ -181,7 +182,6 @@ void PairGranHooke::compute_force(int eflag, int vflag,int addflag)
if (cohesionflag) {
addCohesionForce(i,j,r,Fn_coh);
ccel-=Fn_coh*rinv;
-
}
// relative velocities
diff --git a/src/pair_gran_hooke_history.cpp b/src/pair_gran_hooke_history.cpp
index 3b44172ef4dc240548fef9a185351f4ebdbca982..f5463dd1c1cbb0574a3e28d37a87f69359a3ef40 100644
--- a/src/pair_gran_hooke_history.cpp
+++ b/src/pair_gran_hooke_history.cpp
@@ -675,7 +675,8 @@ void PairGranHookeHistory::init_granular()
if(cohesionflag)
cohEnergyDens1=static_cast<FixPropertyGlobal*>(modify->find_fix_property("cohesionEnergyDensity","property/global","peratomtypepair",max_type,max_type,force->pair_style));
- if(charVelflag) charVel1=static_cast<FixPropertyGlobal*>(modify->find_fix_property("characteristicVelocity","property/global","scalar",0,0,force->pair_style));
+ if(charVelflag)
+ charVel1=static_cast<FixPropertyGlobal*>(modify->find_fix_property("characteristicVelocity","property/global","scalar",0,0,force->pair_style));
//pre-calculate parameters for possible contact material combinations
for(int i=1;i< max_type+1; i++)
@@ -739,7 +740,15 @@ void PairGranHookeHistory::init_granular()
}
}
- if(charVelflag) charVel = charVel1->compute_scalar();
+ if(charVelflag)
+ {
+ charVel = charVel1->compute_scalar();
+ if(sanity_checks)
+ {
+ if(strcmp(update->unit_style,"si") == 0 && charVel < 1e-2)
+ error->all(FLERR,"characteristicVelocity >= 1e-2 required for SI units");
+ }
+ }
// error checks on coarsegraining
if((rollingflag || cohesionflag) && force->cg_active())
diff --git a/src/pointers.h b/src/pointers.h
index 981b5ccd83e710c2151a8b86c52abc80eb34c242..6fce2ed68688048cd3585a65b615ce81170a4052 100644
--- a/src/pointers.h
+++ b/src/pointers.h
@@ -33,6 +33,7 @@ namespace LAMMPS_NS {
#define MIN(A,B) ((A) < (B) ? (A) : (B))
#define MAX(A,B) ((A) > (B) ? (A) : (B))
+#define UNUSED(P) (void)P
class Pointers {
public:
diff --git a/src/region_wedge.cpp b/src/region_wedge.cpp
index 88ab74299861ee74d49df1d82721ff4c0be1da7c..0c7e0af870d7d2faa144ca76e3730ac730a2020f 100644
--- a/src/region_wedge.cpp
+++ b/src/region_wedge.cpp
@@ -193,6 +193,9 @@ RegWedge::RegWedge(LAMMPS *lmp, int narg, char **arg) : Region(lmp, narg, arg)
// of pi_half are incorporated in the bounding box, in case the wedge does
// not reside in one quadrant only.
+ // a ... horizontal direction from viewpoint of axis (cos part of angle)
+ // b ... vertical direction from viewpoint of axis (sin part of angle)
+
double amin, amax;
double bmin, bmax;
@@ -201,29 +204,31 @@ RegWedge::RegWedge(LAMMPS *lmp, int narg, char **arg) : Region(lmp, narg, arg)
// find min and max in both directions for end- and start-angle
amin = MIN(amin, cosang1);
+ amin = MIN(amin, cosang2);
amax = MAX(amax, cosang1);
- amin = MIN(amin, cosang2);
amax = MAX(amax, cosang2);
bmin = MIN(bmin, sinang1);
+ bmin = MIN(bmin, sinang2);
bmax = MAX(bmax, sinang1);
- bmin = MIN(bmin, sinang2);
bmax = MAX(bmax, sinang2);
- // take care of the angles inbetween (suppose angle2-angle1 == pi) then
- // we need to think about what happens in the middle w.r.t. the bbox
-
- double phi = angle1;
- double sinphi, cosphi;
- int n = static_cast<int>(phi/pi_half); // get current multiple of pi_half
-
- while (phi < angle2) {
-
- // round phi to next multiple of pi_half
- n += 1;
+ // sin and cos functions are only monotonous within the four quadrants.
+ // if the wedge crossees quadrands, these functions might have a maximum
+ // there, which should be a limit of the bounding box.
+
+ double phi, sinphi, cosphi;
+ int n = 0;
+ // could be, that angle1 is -190 degrees
+ while (static_cast<double>(n)*pi_half > angle1)
+ n--;
+
+ // since the wedge can have a maximum of 180 degrees
+ // and we start smaller then angle1 i=0 to i=2 suffices
+ for (int i = 0; i < 3; i++, n += pi_half){
phi = static_cast<double>(n)*pi_half;
- // update mins and maxes for this phi
+ if (angle1 <= phi && phi <= angle1 + dang){
sinphi = sin(phi);
cosphi = cos(phi);
@@ -231,10 +236,10 @@ RegWedge::RegWedge(LAMMPS *lmp, int narg, char **arg) : Region(lmp, narg, arg)
amax = MAX(amax, cosphi);
bmin = MIN(bmin, sinphi);
bmax = MAX(bmax, sinphi);
-
+ }
}
- // adjust to incorporate radius
+ // adjust for radius
amin *= radius;
amax *= radius;
bmin *= radius;
diff --git a/src/region_wedge.h b/src/region_wedge.h
index 31644225b34d4e319aa20ca8d869372d67900794..0d0749bb4bd60f65f5dd58a0a889341c14f1e83c 100644
--- a/src/region_wedge.h
+++ b/src/region_wedge.h
@@ -83,7 +83,7 @@ class RegWedge : public Region
double onedivr;
// normal vectors
- // the normal vectors point outside the wedge and are normed
+ // the normal vectors point outside the wedge and are normalized
double normal1[2]; // normal vector to 1st plane of section (angle1)
double normal2[2]; // normal vector to 2nd plane of section (angle2)
diff --git a/src/set.cpp b/src/set.cpp
index 9bcef1f77652f6cbb7be8381cb5965584fd006d3..2001070a32cf11defa846a2ca1ba3e34d4e9f886 100644
--- a/src/set.cpp
+++ b/src/set.cpp
@@ -46,6 +46,7 @@
#include "fix_property_atom.h"
#include "sph_kernels.h"
#include "fix_sph.h"
+#include "update.h"
using namespace LAMMPS_NS;
using namespace MathConst;
@@ -82,6 +83,8 @@ void Set::command(int narg, char **arg)
id = new char[n];
strcpy(id,arg[1]);
select = NULL;
+ add = 0;
+ until = 1;
// loop over keyword/value pairs
// call appropriate routine to reset attributes
@@ -375,6 +378,22 @@ void Set::command(int narg, char **arg)
error->all(FLERR,"Cannot set meso_rho for this atom style");
set(MESO_RHO);
iarg += 2;
+ } else if (strcmp(arg[iarg],"add") == 0){
+ if (iarg+1 > narg)
+ error->all(FLERR,"Illegal set command for add");
+ if(strcmp(arg[iarg+1],"yes") == 0)
+ add = 1;
+ else if(strcmp(arg[iarg+1],"no") == 0)
+ add = 0;
+ else error->all(FLERR,"Illegal 'add' option called");
+ iarg +=2;
+ } else if (strcmp(arg[iarg],"until") == 0){
+ if (iarg+1 > narg)
+ error->all(FLERR,"Illegal set command for until");
+ until = atof(arg[iarg+1]);
+ if (until <= 0.0)
+ error->all(FLERR,"Illegal 'until' option called. Please set keyword value >0");
+ iarg +=2;
} else if (strncmp(arg[iarg],"property/atom",13) == 0) {
if (iarg+1 > narg)
error->all(FLERR,"Illegal set command for property/atom");
@@ -561,17 +580,29 @@ void Set::set(int keyword)
else if (keyword == MESO_RHO) atom->rho[i] = dvalue;
// set desired per-atom property
- else if (keyword == PROPERTYPERATOM) {
+ else if (keyword == PROPERTYPERATOM) {
// if fix was just created, its default values have not been set,
// so have to add a run 0 to call setup
if(updFix->just_created)
error->all(FLERR,"May not use the set command right after fix property/atom without a prior run. Add a 'run 0' between fix property/atom and set");
+ if (add == 0)
+ {
if (updFix->data_style) for (int m = 0; m < nUpdValues; m++)
updFix->array_atom[i][m] = updValues[m];
else updFix->vector_atom[i]=updValues[0];
-
+ }
+ else
+ {
+ currentTimestep = update->ntimestep;
+ if (currentTimestep < until)
+ {
+ if (updFix->data_style) for (int m = 0; m < nUpdValues; m++)
+ updFix->array_atom[i][m] = updValues[m];
+ else updFix->vector_atom[i]=updValues[0];
+ }
+ }
}
// set shape of ellipsoidal particle
diff --git a/src/set.h b/src/set.h
index de8579e4a348db4e0f52a15eb27e49508da77d32..23583698d271f48231f656e6863f1a346e24daa7 100644
--- a/src/set.h
+++ b/src/set.h
@@ -49,11 +49,15 @@ class Set : protected Pointers {
class FixPropertyAtom* updFix;
int nUpdValues;
double *updValues;
+ int add;
+ bigint until;
+ bigint currentTimestep;
void selection(int);
void set(int);
void setrandom(int);
void topology(int);
+
};
}
diff --git a/src/surface_mesh_I.h b/src/surface_mesh_I.h
index a1752686bfc459519eb2ff359e2fb24fe729cc95..7085f6c07c216380ccc02b3eb2bdaac37970a59e 100644
--- a/src/surface_mesh_I.h
+++ b/src/surface_mesh_I.h
@@ -399,7 +399,6 @@ void SurfaceMesh<NUM_NODES,NUM_NEIGH_MAX>::buildNeighbours()
{
int iNode(0), jNode(0), iEdge(0), jEdge(0);
- if(!this->shareNode(i,j,iNode,jNode)) continue;
if(shareEdge(i,j,iEdge,jEdge))
handleSharedEdge(i,iEdge,j,jEdge, areCoplanar(TrackingMesh<NUM_NODES>::id(i),TrackingMesh<NUM_NODES>::id(j)));
@@ -442,7 +441,7 @@ void SurfaceMesh<NUM_NODES,NUM_NEIGH_MAX>::qualityCheck()
int nall = this->sizeLocal()+this->sizeGhost();
int me = this->comm->me;
- // check duplicate elements, n^2 operation
+ // check duplicate elements, n^2/2 operation
for(int i = 0; i < nlocal; i++)
{
@@ -490,7 +489,8 @@ void SurfaceMesh<NUM_NODES,NUM_NEIGH_MAX>::qualityCheck()
fprintf(this->screen,"Mesh %s: %d mesh elements have more than %d neighbors \n",
this->mesh_id_,this->nTooManyNeighs(),NUM_NEIGH_MAX);
- this->error->one(FLERR,"Fix mesh: Bad mesh, cannot continue. Possibly corrupt elements with too many neighbors");
+ this->error->one(FLERR,"Fix mesh: Bad mesh, cannot continue. Possibly corrupt elements with too many neighbors.\n"
+ "If you know what you're doing, you can try to change the definition of SurfaceMeshBase in tri_mesh.h and recompile");
}
if(nOverlapping() > 0)
@@ -578,7 +578,7 @@ bool SurfaceMesh<NUM_NODES,NUM_NEIGH_MAX>::areCoplanar(int tag_a, int tag_b)
double dot = vectorDot3D(surfaceNorm(a),surfaceNorm(b));
// need fabs in case surface normal is other direction
- if(fabs(dot) > curvature_) return true;
+ if(fabs(dot) >= curvature_) return true;
else return false;
}
@@ -714,23 +714,26 @@ void SurfaceMesh<NUM_NODES,NUM_NEIGH_MAX>::growSurface(int iSrf, double by)
template<int NUM_NODES, int NUM_NEIGH_MAX>
bool SurfaceMesh<NUM_NODES,NUM_NEIGH_MAX>::shareEdge(int iSrf, int jSrf, int &iEdge, int &jEdge)
{
- int i,j;
- if(this->shareNode(iSrf,jSrf,i,j)){
+ int iNode1,jNode1,iNode2,jNode2;
+ if(this->share2Nodes(iSrf,jSrf,iNode1,jNode1,iNode2,jNode2)){
// following implementation of shareNode(), the only remaining option to
- // share an edge is that the next node of iSrf is equal to the previous
+ // share an edge is that the next node of iSrf is equal to the next or previous
// node if jSrf
- if(i==0 && MultiNodeMesh<NUM_NODES>::nodesAreEqual(iSrf,NUM_NODES-1,jSrf,(j+1)%NUM_NODES)){
- iEdge = NUM_NODES-1;
- jEdge = j;
- return true;
- }
- if(MultiNodeMesh<NUM_NODES>::nodesAreEqual
- (iSrf,(i+1)%NUM_NODES,jSrf,(j-1+NUM_NODES)%NUM_NODES)){
- iEdge = i;//(ii-1+NUM_NODES)%NUM_NODES;
- jEdge = (j-1+NUM_NODES)%NUM_NODES;
- return true;
- }
+
+ if(2 == iNode1+iNode2)
+ iEdge = 2;
+
+ else
+ iEdge = MathExtraLiggghts::min(iNode1,iNode2);
+
+ if(2 == jNode1+jNode2)
+ jEdge = 2;
+ else
+ jEdge = MathExtraLiggghts::min(jNode1,jNode2);
+
+ return true;
}
+
iEdge = -1; jEdge = -1;
return false;
}
@@ -741,6 +744,7 @@ template<int NUM_NODES, int NUM_NEIGH_MAX>
void SurfaceMesh<NUM_NODES,NUM_NEIGH_MAX>::handleSharedEdge(int iSrf, int iEdge, int jSrf, int jEdge,
bool coplanar, bool neighflag)
{
+
if(neighflag)
{
if(nNeighs_(iSrf) == NUM_NEIGH_MAX || nNeighs_(jSrf) == NUM_NEIGH_MAX)
@@ -784,7 +788,7 @@ void SurfaceMesh<NUM_NODES,NUM_NEIGH_MAX>::handleSharedEdge(int iSrf, int iEdge,
edgeActive(jSrf)[jEdge] = false;
}
}
- else
+ else // coplanar
{
if(!coplanar) this->error->one(FLERR,"internal error");
@@ -828,10 +832,9 @@ int SurfaceMesh<NUM_NODES,NUM_NEIGH_MAX>::handleCorner(int iSrf, int iNode,
}
}
- // deactivate all
if(hasTwoColinearEdges || !anyActiveEdge)
cornerActive(iSrf)[iNode] = false;
- // let the highest ID live
+
else if(TrackingMesh<NUM_NODES>::id(iSrf) == maxId)
cornerActive(iSrf)[iNode] = true;
else
diff --git a/src/vector_liggghts.h b/src/vector_liggghts.h
index a8f17b54fc5d6013207ce298f0ba6745ee771e96..2ca74f68a6ce8b892b0651d5807a0062236612a5 100644
--- a/src/vector_liggghts.h
+++ b/src/vector_liggghts.h
@@ -48,7 +48,7 @@ inline void vectorConstruct3D(int *v,int x, int y, int z)
inline void vectorNormalize3D(double *v)
{
double norm = sqrt(v[0]*v[0]+v[1]*v[1]+v[2]*v[2]);
- double invnorm = (norm == 0) ? 0. : 1./norm;
+ double invnorm = (norm == 0.) ? 0. : 1./norm;
v[0] *= invnorm;
v[1] *= invnorm;
v[2] *= invnorm;
@@ -110,6 +110,13 @@ inline void vectorCopy3D(double *from,double *to)
to[2]=from[2];
}
+inline void vectorFlip3D(double *v)
+{
+ v[0]=-v[0];
+ v[1]=-v[1];
+ v[2]=-v[2];
+}
+
inline void vectorCopyN(double *from,double *to,int N)
{
for(int i = 0; i < N; i++)
@@ -281,6 +288,12 @@ inline void vectorSubtract3D(const double *v1,const double *v2, double *result)
result[2]=v1[2]-v2[2];
}
+inline void vectorSubtract2D(const double *v1,const double *v2, double *result)
+{
+ result[0]=v1[0]-v2[0];
+ result[1]=v1[1]-v2[1];
+}
+
inline void vectorCross3D(const double *v1,const double *v2, double *result)
{
result[0]=v1[1]*v2[2]-v1[2]*v2[1];
@@ -288,6 +301,15 @@ inline void vectorCross3D(const double *v1,const double *v2, double *result)
result[2]=v1[0]*v2[1]-v1[1]*v2[0];
}
+inline double vectorCrossMag3D(const double *v1,const double *v2)
+{
+ double res[3];
+ res[0]=v1[1]*v2[2]-v1[2]*v2[1];
+ res[1]=v1[2]*v2[0]-v1[0]*v2[2];
+ res[2]=v1[0]*v2[1]-v1[1]*v2[0];
+ return vectorMag3D(res);
+}
+
inline void vectorZeroize3D(double *v)
{
v[0]=0.;
@@ -400,6 +422,11 @@ inline void bufToVector4D(double *vec,double *buf,int &m)
vec[3] = buf[m++];
}
+inline void printVec2D(FILE *out, const char *name, double *vec)
+{
+ fprintf(out," vector %s: %e %e\n",name,vec[0],vec[1]);
+}
+
inline void printVec3D(FILE *out, const char *name, double *vec)
{
fprintf(out," vector %s: %e %e %e\n",name,vec[0],vec[1],vec[2]);
diff --git a/src/verlet_implicit.cpp b/src/verlet_implicit.cpp
index 758600f9abac35b482a4ef17d87062e1ce496cff..dbaa011ead8a6a8fb66a0b53130bc2e84ded9eb1 100644
--- a/src/verlet_implicit.cpp
+++ b/src/verlet_implicit.cpp
@@ -73,25 +73,11 @@ void VerletImplicit::run(int n)
ntimestep = ++update->ntimestep;
- do
- {
- ev_set(ntimestep);
-
- // initial time integration
-
- modify->initial_integrate(vflag);
-
- if (n_post_integrate) modify->post_integrate();
+ // neigh list build if required
- // regular communication vs neighbor list rebuild
+ nflag = neighbor->decide();
- nflag = neighbor->decide();
-
- if (nflag == 0) {
- timer->stamp();
- comm->forward_comm();
- timer->stamp(TIME_COMM);
- } else {
+ if (nflag == 1) {
if (n_pre_exchange) modify->pre_exchange();
@@ -116,7 +102,23 @@ void VerletImplicit::run(int n)
neighbor->build();
timer->stamp(TIME_NEIGHBOR);
- }
+ }
+
+ do
+ {
+ ev_set(ntimestep);
+
+ // initial time integration
+
+ modify->initial_integrate(vflag);
+
+ if (n_post_integrate) modify->post_integrate();
+
+ // regular communication here
+
+ timer->stamp();
+ comm->forward_comm();
+ timer->stamp(TIME_COMM);
// force computations
@@ -155,7 +157,6 @@ void VerletImplicit::run(int n)
if (n_post_force) modify->post_force(vflag);
modify->final_integrate();
-
}
while(modify->iterate_implicitly());
diff --git a/src/version_liggghts.h b/src/version_liggghts.h
index 86e6ee1362609edc7d349ac857e8510e682591a7..53eb557bbf13d28d8facf7acebcea4ca88d46923 100644
--- a/src/version_liggghts.h
+++ b/src/version_liggghts.h
@@ -1 +1 @@
-#define LIGGGHTS_VERSION "LIGGGHTS-PUBLIC 2.3.7, compiled 2013-09-06-15:08:36 by ckloss"
+#define LIGGGHTS_VERSION "LIGGGHTS-PUBLIC 2.3.8, compiled 2013-10-10-16:13:21 by ckloss"
diff --git a/src/version_liggghts.txt b/src/version_liggghts.txt
index 00355e29d11ac4a7ee62410face3b0adb50201ac..bc4abe86dedf538f8c5783cf2cb2b7dcb5188deb 100644
--- a/src/version_liggghts.txt
+++ b/src/version_liggghts.txt
@@ -1 +1 @@
-2.3.7
+2.3.8
diff --git a/src/volume_mesh_I.h b/src/volume_mesh_I.h
index fe4bafd027ed0eaed33c29bb193104cf81e4d3c6..20fae8b3e4d0010a7a7b663966a140161236966c 100644
--- a/src/volume_mesh_I.h
+++ b/src/volume_mesh_I.h
@@ -346,7 +346,7 @@ void VolumeMesh<NUM_NODES,NUM_FACES,NUM_NODES_PER_FACE>::buildNeighbours()
{
int iNode(0), jNode(0), iEdge(0), jEdge(0);
- if(!this->shareNode(i,j,iNode,jNode)) continue;
+ if(0 == this->nSharedNodes(i,j)) continue;
if(shareFace(i,j,iFace,jFace))
{