[liggghts] 01/03: Imported Upstream version 2.3.8
Anton Gladky
gladk at moszumanska.debian.org
Wed Nov 20 15:28:38 UTC 2013
This is an automated email from the git hooks/post-receive script.
gladk pushed a commit to branch master
in repository liggghts.
commit 3330f02ddffab3c8dcb8cedbbe9587400fe587a9
Author: Anton Gladky <gladk at debian.org>
Date: Sun Nov 10 20:22:52 2013 +0100
Imported Upstream version 2.3.8
---
doc/Eqs/conduction_rate.png | Bin 0 -> 3548 bytes
doc/Eqs/model_A.png | Bin 0 -> 4328 bytes
doc/Eqs/model_C_part_1.png | Bin 0 -> 12830 bytes
doc/Eqs/model_C_part_2.png | Bin 0 -> 5540 bytes
doc/Eqs/model_C_part_3.png | Bin 0 -> 9426 bytes
doc/Eqs/model_b1.png | Bin 0 -> 28056 bytes
doc/Eqs/model_b2_part1.png | Bin 0 -> 4992 bytes
doc/Eqs/model_b2_part2.png | Bin 0 -> 5999 bytes
doc/Eqs/modelb1.png | Bin 0 -> 5445 bytes
doc/Eqs/scsr.png | Bin 0 -> 6372 bytes
doc/Eqs/variables.png | Bin 0 -> 38743 bytes
doc/Eqs/vb.png | Bin 0 -> 3087 bytes
doc/Eqs/vbi.png | Bin 0 -> 7242 bytes
doc/Eqs/vbj.png | Bin 0 -> 7289 bytes
doc/Manual.pdf | Bin 6550328 -> 6551572 bytes
doc/fix_massflow_mesh.html | 16 +-
doc/fix_massflow_mesh.txt | 16 +-
doc/fix_property.html | 41 +++-
doc/fix_property.txt | 40 +++-
doc/neigh_modify.html | 8 +-
doc/neigh_modify.txt | 8 +-
doc/set.html | 10 +
doc/set.txt | 10 +
.../Tutorials_public/heatTransfer_1/in.heatGran | 2 +-
.../Tutorials_public/heatTransfer_1/log.liggghts | 202 --------------------
.../Tutorials_public/heatTransfer_2/in.heatGran | 2 +-
src/atom_vec_ellipsoid.cpp | 1 +
src/atom_vec_sphere.cpp | 14 ++
src/atom_vec_sphere.h | 2 +
...{domain_wedge_dummy.h => atom_vec_sphere_w.cpp} | 31 ++-
src/comm.cpp | 103 ++++++++--
src/comm.h | 8 +
src/comm_I.h | 81 ++++++++
src/compute_erotate_sphere.cpp | 7 +-
src/compute_erotate_sphere.h | 1 +
src/compute_ke.cpp | 6 +-
src/compute_ke.h | 1 +
src/compute_ke_atom.cpp | 5 +-
src/compute_ke_atom.h | 1 +
src/compute_reduce.cpp | 20 +-
src/container_base.cpp | 16 +-
src/container_base.h | 7 +-
src/container_base_I.h | 22 +--
src/debug_liggghts.h | 1 +
src/domain.cpp | 39 +---
src/domain.h | 115 ++---------
src/domain_I.h | 178 +++++++++++++++++
src/domain_wedge_dummy.h | 14 ++
src/dump_mesh_vtk.cpp | 6 +-
src/fix.h | 5 +-
src/fix_ave_spatial.cpp | 6 +-
src/fix_cfd_coupling_force.cpp | 2 +-
src/fix_cfd_coupling_force_implicit.cpp | 126 +++++++++++-
src/fix_cfd_coupling_force_implicit.h | 5 +
src/fix_contact_history.cpp | 43 ++---
src/fix_contact_history.h | 6 +
src/fix_heat_gran.cpp | 1 -
src/fix_insert.cpp | 8 +-
src/fix_insert_pack.cpp | 3 +-
src/fix_massflow_mesh.cpp | 113 +++++++++--
src/fix_massflow_mesh.h | 10 +
src/fix_mesh.cpp | 4 +-
src/fix_mesh_surface_stress.cpp | 37 ++--
src/fix_mesh_surface_stress.h | 1 +
src/fix_mesh_surface_stress_servo.cpp | 16 +-
src/fix_move.cpp | 6 +
src/fix_neighlist_mesh.cpp | 2 +
src/fix_property_atom_tracer.cpp | 8 +-
src/fix_property_atom_tracer.h | 3 +-
src/fix_property_global.cpp | 1 +
src/fix_scalar_transport_equation.cpp | 2 +-
src/general_container.h | 2 +-
src/general_container_I.h | 7 +-
src/lammps.cpp | 3 +-
src/math_extra_liggghts.h | 8 +-
src/modify.cpp | 5 +
src/modify.h | 1 +
src/modify_liggghts.cpp | 16 +-
src/multi_node_mesh.h | 8 +-
src/multi_node_mesh_I.h | 34 +++-
src/multi_node_mesh_parallel_buffer_I.h | 30 +--
src/multi_vector_container.h | 6 +-
src/neigh_gran.cpp | 12 +-
src/neighbor.cpp | 12 +-
src/pair_gran.cpp | 3 +
src/pair_gran_hooke.cpp | 2 +-
src/pair_gran_hooke_history.cpp | 17 +-
src/pair_hybrid.cpp | 23 +++
src/pair_hybrid.h | 1 +
src/pair_sph.cpp | 22 +++
src/pair_sph.h | 3 +
src/pair_sph_artvisc_tenscorr.cpp | 4 +
src/pointers.h | 1 +
src/primitive_wall.h | 2 +-
src/read_data.cpp | 90 +++++----
src/region.cpp | 7 +-
src/region_wedge.cpp | 37 ++--
src/region_wedge.h | 2 +-
src/scalar_container.h | 6 +-
src/set.cpp | 35 +++-
src/set.h | 4 +
src/surface_mesh_I.h | 45 +++--
src/vector_container.h | 6 +-
src/vector_liggghts.h | 29 ++-
src/verlet_implicit.cpp | 39 ++--
src/version_liggghts.h | 2 +-
src/version_liggghts.txt | 2 +-
src/volume_mesh_I.h | 2 +-
108 files changed, 1288 insertions(+), 682 deletions(-)
diff --git a/doc/Eqs/conduction_rate.png b/doc/Eqs/conduction_rate.png
new file mode 100644
index 0000000..391ea8c
Binary files /dev/null and b/doc/Eqs/conduction_rate.png differ
diff --git a/doc/Eqs/model_A.png b/doc/Eqs/model_A.png
new file mode 100644
index 0000000..f84de96
Binary files /dev/null and b/doc/Eqs/model_A.png differ
diff --git a/doc/Eqs/model_C_part_1.png b/doc/Eqs/model_C_part_1.png
new file mode 100644
index 0000000..ec6e131
Binary files /dev/null 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
new file mode 100644
index 0000000..e8030e8
Binary files /dev/null and b/doc/Eqs/model_C_part_2.png differ
diff --git a/doc/Eqs/model_C_part_3.png b/doc/Eqs/model_C_part_3.png
new file mode 100644
index 0000000..e79e2a6
Binary files /dev/null and b/doc/Eqs/model_C_part_3.png differ
diff --git a/doc/Eqs/model_b1.png b/doc/Eqs/model_b1.png
new file mode 100644
index 0000000..c993913
Binary files /dev/null and b/doc/Eqs/model_b1.png differ
diff --git a/doc/Eqs/model_b2_part1.png b/doc/Eqs/model_b2_part1.png
new file mode 100644
index 0000000..c072a9c
Binary files /dev/null and b/doc/Eqs/model_b2_part1.png differ
diff --git a/doc/Eqs/model_b2_part2.png b/doc/Eqs/model_b2_part2.png
new file mode 100644
index 0000000..91eb05a
Binary files /dev/null and b/doc/Eqs/model_b2_part2.png differ
diff --git a/doc/Eqs/modelb1.png b/doc/Eqs/modelb1.png
new file mode 100644
index 0000000..0540687
Binary files /dev/null and b/doc/Eqs/modelb1.png differ
diff --git a/doc/Eqs/scsr.png b/doc/Eqs/scsr.png
new file mode 100644
index 0000000..3367504
Binary files /dev/null and b/doc/Eqs/scsr.png differ
diff --git a/doc/Eqs/variables.png b/doc/Eqs/variables.png
new file mode 100644
index 0000000..e389d77
Binary files /dev/null and b/doc/Eqs/variables.png differ
diff --git a/doc/Eqs/vb.png b/doc/Eqs/vb.png
new file mode 100644
index 0000000..50ea223
Binary files /dev/null and b/doc/Eqs/vb.png differ
diff --git a/doc/Eqs/vbi.png b/doc/Eqs/vbi.png
new file mode 100644
index 0000000..f08f338
Binary files /dev/null and b/doc/Eqs/vbi.png differ
diff --git a/doc/Eqs/vbj.png b/doc/Eqs/vbj.png
new file mode 100644
index 0000000..6d461a5
Binary files /dev/null and b/doc/Eqs/vbj.png differ
diff --git a/doc/Manual.pdf b/doc/Manual.pdf
index 40e674c..f86b722 100644
Binary files a/doc/Manual.pdf and b/doc/Manual.pdf differ
diff --git a/doc/fix_massflow_mesh.html b/doc/fix_massflow_mesh.html
index 3dd47e4..a87ab77 100644
--- a/doc/fix_massflow_mesh.html
+++ b/doc/fix_massflow_mesh.html
@@ -31,7 +31,7 @@
<LI>zero or more keyword/value pairs may be appended to args
-<LI>keywords = <I>count</I> or <I>point_at_outlet</I> or <I>append</I> or <I>file</I> or <I>screen</I>
+<LI>keywords = <I>count</I> or <I>point_at_outlet</I> or <I>append</I> or <I>file</I> or <I>screen</I> or <I>delete_atoms</I>
<PRE> <I>count</I> value = <I>once</I> or <I>multiple</I>
once = count particles only once
@@ -43,7 +43,9 @@
<I>file</I> value = filename
<I>append</I> value = filename
filename = name of the file to print radius, position and velocity values of the particles
- <I>screen</I> value = <I>yes</I> or <I>no</I>
+ <I>screen</I> value = <I>yes</I> or <I>no</I>
+ <I>delete_atoms</I> value = <I>yes</I> or <I>no</I>
+ yes = to remove the particles that pass through the mesh surface
</PRE>
</UL>
@@ -94,7 +96,8 @@ by specifying a filename.
<P>If the <I>screen</I> keyword is used, output by this fix to the screen and
logfile can be turned on or off as desired.
</P>
-<P>If the <I>delete_atoms</I> keyword is used then the particles passing through the mesh surface are deleted.
+<P>If the <I>delete_atoms</I> keyword is used then the particles passing through the mesh
+surface are deleted at the next re-neighboring step.
</P>
<P><B>Restart, fix_modify, output, run start/stop, minimize info:</B>
</P>
@@ -103,14 +106,15 @@ logfile can be turned on or off as desired.
<P>This fix computes a per-atom vector (the marker) which can be accessed by
various <A HREF = "Section_howto.html#howto_15">output commands</A>. The per-atom vector
(i.e., the marker) can be accessed by dumps by using "f_massflow_ID", .
-This fix also computes a global vector of length 4. This vector can be
+This fix also computes a global vector of length 6. This vector can be
accessed via "f_ID", where ID is the fix id. The first vector component
is equal to the total mass which has crossed the mesh surface, the second vector
component indicates the particle count. The third vector component
is equal to the total mass which has crossed the mesh surface since the last output
divived by the time since the last output (i.e., the mass flow rate), the fourth vector
component indicates the particle count since the last output divived by the time
-since the last output (i.e., the number rate of particles). This vector
+since the last output (i.e., the number rate of particles). The fifth and sixth vector
+components are the deleted mass and the number of deleted particles. This vector
can also be accessed by various <A HREF = "Section_howto.html#howto_15">output commands</A>.
</P>
<P><B>Restrictions:</B>
@@ -123,6 +127,6 @@ can also be accessed by various <A HREF = "Section_howto.html#howto_15">output c
</P>
<P><B>Default:</B>
</P>
-<P><I>count</I> = multiple
+<P><I>count</I> = multiple, <I>inside_out</I> =false, <I>delete_atoms</I> = false
</P>
</HTML>
diff --git a/doc/fix_massflow_mesh.txt b/doc/fix_massflow_mesh.txt
index 115bc5f..60cc18e 100644
--- a/doc/fix_massflow_mesh.txt
+++ b/doc/fix_massflow_mesh.txt
@@ -20,7 +20,7 @@ mesh-ID = ID of a "fix mesh/surface"_fix_mesh_surface.html command :l
vec_side = obligatory keyword :l
vx, vy, vz = vector components defining the "outside" of the mesh :l
zero or more keyword/value pairs may be appended to args :l
-keywords = {count} or {point_at_outlet} or {append} or {file} or {screen} :l
+keywords = {count} or {point_at_outlet} or {append} or {file} or {screen} or {delete_atoms} :l
{count} value = {once} or {multiple}
once = count particles only once
multiple = allow particles to be counted multiple times
@@ -31,7 +31,9 @@ keywords = {count} or {point_at_outlet} or {append} or {file} or {screen} :l
{file} value = filename
{append} value = filename
filename = name of the file to print radius, position and velocity values of the particles
- {screen} value = {yes} or {no} :pre
+ {screen} value = {yes} or {no}
+ {delete_atoms} value = {yes} or {no}
+ yes = to remove the particles that pass through the mesh surface :pre
:ule
@@ -81,7 +83,8 @@ by specifying a filename.
If the {screen} keyword is used, output by this fix to the screen and
logfile can be turned on or off as desired.
-If the {delete_atoms} keyword is used then the particles passing through the mesh surface are deleted.
+If the {delete_atoms} keyword is used then the particles passing through the mesh
+surface are deleted at the next re-neighboring step.
[Restart, fix_modify, output, run start/stop, minimize info:]
@@ -90,14 +93,15 @@ Information about this fix is written to "binary restart files"_restart.html .
This fix computes a per-atom vector (the marker) which can be accessed by
various "output commands"_Section_howto.html#howto_15. The per-atom vector
(i.e., the marker) can be accessed by dumps by using "f_massflow_ID", .
-This fix also computes a global vector of length 4. This vector can be
+This fix also computes a global vector of length 6. This vector can be
accessed via "f_ID", where ID is the fix id. The first vector component
is equal to the total mass which has crossed the mesh surface, the second vector
component indicates the particle count. The third vector component
is equal to the total mass which has crossed the mesh surface since the last output
divived by the time since the last output (i.e., the mass flow rate), the fourth vector
component indicates the particle count since the last output divived by the time
-since the last output (i.e., the number rate of particles). This vector
+since the last output (i.e., the number rate of particles). The fifth and sixth vector
+components are the deleted mass and the number of deleted particles. This vector
can also be accessed by various "output commands"_Section_howto.html#howto_15.
[Restrictions:]
@@ -110,4 +114,4 @@ Currently, this feature does not support multi-sphere particles.
[Default:]
-{count} = multiple
+{count} = multiple, {inside_out} =false, {delete_atoms} = false
diff --git a/doc/fix_property.html b/doc/fix_property.html
index 01e2c18..ad43085 100644
--- a/doc/fix_property.html
+++ b/doc/fix_property.html
@@ -55,17 +55,40 @@ fix id group property/global variablename style stylearg defaultvalue(s)...
</UL>
+<P><B>Examples:</B>
+</P>
+<PRE>fix m3 all property/global coefficientRestitution peratomtypepair 1 0.3
+fix m5 all property/global characteristicVelocity scalar 2.
+fix uf all property/atom uf vector yes no no 0. 0. 0.
+</PRE>
<P><B>LIGGGHTS vs. LAMMPS Info:</B>
</P>
<P>This LIGGGHTS command is not available in LAMMPS.
</P>
<P><B>Description:</B>
</P>
-<P><B>Fix property/atom</B> reserves per-atom properties to be accessed by the user or other fixes. Style scalar reserves one value per atom, style vector multiple values per atoms, where the number of defaultvalues (that are assigned to the atoms at creation) determines the length of the vector . The group of atoms the fix is applied to is always "all", irrespective of which group is used for the fix command . If you want to assign different values for different groups, you can use the [...]
-</P>
-<P><B>Fix property/global</B> reserves global properties to be accessed by the user or other fixes or pair styles. The number of defaultvalues determines the length of the vector / the number of matrix components . For style vector or atomtype, the user provides the number of vector components . For style matrix or atomtypepair, the user provides the number of matrix columns (nCols) .
-</P>
-<P>Example: nCols= 2 and defaultvalues = 1 2 3 4 5 6 would be mapped into a matrix like
+<P><B>Fix property/atom</B> reserves per-atom properties to be accessed by the user or other fixes.
+Style <I>scalar</I> reserves one value per atom, style <I>vector</I> multiple values per atoms, where
+the number of <I>defaultvalues</I> (that are assigned to the atoms at creation) determines the
+length of the vector. The group of atoms the fix is applied to is always "all", irrespective
+of which group is used for the fix command . If you want to assign different values for
+different groups, you can use the <A HREF = "set.html">set</A> command with keyword 'property/atom'.
+Keyword <I>restartvalues</I> determines whether information about the values stored by this fix
+is written to binary restart files.
+Keyword <I>communicate_ghost_value</I> determines whether information about the values stored by this fix
+can be communicated to ghost particles (forward communication). The exact location during a time-step
+when this happens depends on the model that uses this fix.
+Keyword <I>communicate_reverse_ghost_value</I> determines whether information about the values stored by this fix
+can be communicated from ghost particles to owned particles (reverse communication). The exact location
+during a time-step when this happens depends on the model that uses this fix.
+</P>
+<P><B>Fix property/global</B> reserves global properties to be accessed by the user or other
+fixes or pair styles. The number of defaultvalues determines the length of the vector /
+the number of matrix components . For style <I>vector</I> or <I>atomtype</I>, the user provides
+the number of vector components . For style <I>matrix</I> or <I>atomtypepair</I>, the user provides
+the number of matrix columns (<I>nCols</I>) .
+</P>
+<P>Example: <I>nCols</I>= 2 and <I>defaultvalues</I> = 1 2 3 4 5 6 would be mapped into a matrix like
</P>
<P>[ 1 2 ]
</P>
@@ -73,15 +96,15 @@ fix id group property/global variablename style stylearg defaultvalue(s)...
</P>
<P>[ 5 6 ]
</P>
-<P>Note that the number of default values must thus be a multiple of nCols.
-Note that vector and atomtype do the same thing, atomtype is just provided to make input scripts more readable .
-Note that matrix and atomtypepair both refer to a matrix of global values. However, a matrix defined via 'atomtypepair' is required to be symmetric.
+<P>Note that the number of default values must thus be a multiple of <I>nCols</I>.
+Note that <I>vector</I> and <I>atomtype</I> do the same thing, <I>atomtype</I> is just provided to make input scripts more readable .
+Note that <I>matrix</I> and <I>atomtypepair</I> both refer to a matrix of global values. However, a matrix defined via <I>atomtypepair</I> is required to be symmetric.
</P>
<P>Note that the group of atoms the fix is applied to is ignored (as the fix is not applied to atoms, but defines values of global scope).
</P>
<P><B>Restart, fix_modify, output, run start/stop, minimize info:</B>
</P>
-<P>Information about this fix is written to <A HREF = "restart.html">binary restart files</A> if you set 'restartvalue' to 'yes'.
+<P>Information about this fix is written to <A HREF = "restart.html">binary restart files</A> if you set <I>restartvalue</I> to 'yes'.
</P>
<P><B>Restrictions:</B> none
</P>
diff --git a/doc/fix_property.txt b/doc/fix_property.txt
index 09928fa..df9e9d4 100644
--- a/doc/fix_property.txt
+++ b/doc/fix_property.txt
@@ -38,6 +38,11 @@ fix property/atom:
communicate_reverse_ghost_value = yes or no :l
:ule
+[Examples:]
+
+fix m3 all property/global coefficientRestitution peratomtypepair 1 0.3
+fix m5 all property/global characteristicVelocity scalar 2.
+fix uf all property/atom uf vector yes no no 0. 0. 0. :pre
[LIGGGHTS vs. LAMMPS Info:]
@@ -45,11 +50,28 @@ This LIGGGHTS command is not available in LAMMPS.
[Description:]
-[Fix property/atom] reserves per-atom properties to be accessed by the user or other fixes. Style scalar reserves one value per atom, style vector multiple values per atoms, where the number of defaultvalues (that are assigned to the atoms at creation) determines the length of the vector . The group of atoms the fix is applied to is always "all", irrespective of which group is used for the fix command . If you want to assign different values for different groups, you can use the set comm [...]
-
-[Fix property/global] reserves global properties to be accessed by the user or other fixes or pair styles. The number of defaultvalues determines the length of the vector / the number of matrix components . For style vector or atomtype, the user provides the number of vector components . For style matrix or atomtypepair, the user provides the number of matrix columns (nCols) .
-
-Example: nCols= 2 and defaultvalues = 1 2 3 4 5 6 would be mapped into a matrix like
+[Fix property/atom] reserves per-atom properties to be accessed by the user or other fixes.
+Style {scalar} reserves one value per atom, style {vector} multiple values per atoms, where
+the number of {defaultvalues} (that are assigned to the atoms at creation) determines the
+length of the vector. The group of atoms the fix is applied to is always "all", irrespective
+of which group is used for the fix command . If you want to assign different values for
+different groups, you can use the "set"_set.html command with keyword 'property/atom'.
+Keyword {restartvalues} determines whether information about the values stored by this fix
+is written to binary restart files.
+Keyword {communicate_ghost_value} determines whether information about the values stored by this fix
+can be communicated to ghost particles (forward communication). The exact location during a time-step
+when this happens depends on the model that uses this fix.
+Keyword {communicate_reverse_ghost_value} determines whether information about the values stored by this fix
+can be communicated from ghost particles to owned particles (reverse communication). The exact location
+during a time-step when this happens depends on the model that uses this fix.
+
+[Fix property/global] reserves global properties to be accessed by the user or other
+fixes or pair styles. The number of defaultvalues determines the length of the vector /
+the number of matrix components . For style {vector} or {atomtype}, the user provides
+the number of vector components . For style {matrix} or {atomtypepair}, the user provides
+the number of matrix columns ({nCols}) .
+
+Example: {nCols}= 2 and {defaultvalues} = 1 2 3 4 5 6 would be mapped into a matrix like
\[ 1 2 \]
@@ -57,15 +79,15 @@ Example: nCols= 2 and defaultvalues = 1 2 3 4 5 6 would be mapped into a matrix
\[ 5 6 \]
-Note that the number of default values must thus be a multiple of nCols.
-Note that vector and atomtype do the same thing, atomtype is just provided to make input scripts more readable .
-Note that matrix and atomtypepair both refer to a matrix of global values. However, a matrix defined via 'atomtypepair' is required to be symmetric.
+Note that the number of default values must thus be a multiple of {nCols}.
+Note that {vector} and {atomtype} do the same thing, {atomtype} is just provided to make input scripts more readable .
+Note that {matrix} and {atomtypepair} both refer to a matrix of global values. However, a matrix defined via {atomtypepair} is required to be symmetric.
Note that the group of atoms the fix is applied to is ignored (as the fix is not applied to atoms, but defines values of global scope).
[Restart, fix_modify, output, run start/stop, minimize info:]
-Information about this fix is written to "binary restart files"_restart.html if you set 'restartvalue' to 'yes'.
+Information about this fix is written to "binary restart files"_restart.html if you set {restartvalue} to 'yes'.
[Restrictions:] none
diff --git a/doc/neigh_modify.html b/doc/neigh_modify.html
index b0b128a..3565337 100644
--- a/doc/neigh_modify.html
+++ b/doc/neigh_modify.html
@@ -43,6 +43,8 @@
N = number of pairs stored in a single neighbor page
<I>one</I> value = N
N = max number of neighbors of one atom
+ <I>contactHistoryDistanceFactor</I> value = N
+ N = contact history distance factor used to keep contact history of non-contacting particles (must be > 1).
<I>binsize</I> value = size
size = bin size for neighbor list construction (distance units)
</PRE>
@@ -54,7 +56,8 @@
neigh_modify exclude type 2 3
neigh_modify exclude group frozen frozen check no
neigh_modify exclude group residue1 chain3
-neigh_modify exclude molecule rigid
+neigh_modify exclude molecule rigid
+neigh_modify delay 0 contactHistoryDistanceFactor 2.0
</PRE>
<P><B>Description:</B>
</P>
@@ -80,6 +83,9 @@ lists should be rebuilt.
<P>When the rRESPA integrator is used (see the <A HREF = "run_style.html">run_style</A>
command), the <I>every</I> and <I>delay</I> parameters refer to the longest
(outermost) timestep.
+The <I>contactHistoryDistanceFactor</I> setting has to be used while using fix liquid transfer.
+This factor must be either 1 or greater than 1, to ensure that the extra contact history values are not lost
+when the neighbor lists are rebuilt.
</P>
<P>The <I>include</I> option limits the building of pairwise neighbor lists to
atoms in the specified group. This can be useful for models where a
diff --git a/doc/neigh_modify.txt b/doc/neigh_modify.txt
index 51282b2..0bf8511 100644
--- a/doc/neigh_modify.txt
+++ b/doc/neigh_modify.txt
@@ -39,6 +39,8 @@ keyword = {delay} or {every} or {check} or {once} or {include} or {exclude} or {
N = number of pairs stored in a single neighbor page
{one} value = N
N = max number of neighbors of one atom
+ {contactHistoryDistanceFactor} value = N
+ N = contact history distance factor used to keep contact history of non-contacting particles (must be > 1).
{binsize} value = size
size = bin size for neighbor list construction (distance units) :pre
:ule
@@ -49,7 +51,8 @@ neigh_modify every 2 delay 10 check yes page 100000
neigh_modify exclude type 2 3
neigh_modify exclude group frozen frozen check no
neigh_modify exclude group residue1 chain3
-neigh_modify exclude molecule rigid :pre
+neigh_modify exclude molecule rigid
+neigh_modify delay 0 contactHistoryDistanceFactor 2.0 :pre
[Description:]
@@ -75,6 +78,9 @@ lists should be rebuilt.
When the rRESPA integrator is used (see the "run_style"_run_style.html
command), the {every} and {delay} parameters refer to the longest
(outermost) timestep.
+The {contactHistoryDistanceFactor} setting has to be used while using fix liquid transfer.
+This factor must be either 1 or greater than 1, to ensure that the extra contact history values are not lost
+when the neighbor lists are rebuilt.
The {include} option limits the building of pairwise neighbor lists to
atoms in the specified group. This can be useful for models where a
diff --git a/doc/set.html b/doc/set.html
index 6cfebeb..34a9b49 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 5ab92b7..0bcc0ff 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/examples/LIGGGHTS/Tutorials_public/heatTransfer_1/in.heatGran b/examples/LIGGGHTS/Tutorials_public/heatTransfer_1/in.heatGran
index aad8494..19e9e1f 100644
--- a/examples/LIGGGHTS/Tutorials_public/heatTransfer_1/in.heatGran
+++ b/examples/LIGGGHTS/Tutorials_public/heatTransfer_1/in.heatGran
@@ -49,7 +49,7 @@ fix pts1 all particletemplate/sphere 1 atom_type 1 density constant 8000 radius
fix pdd1 all particledistribution/discrete 1. 1 pts1 1.0
fix ins all insert/pack seed 100001 distributiontemplate pdd1 vel constant 0. 0. -0.3 &
- insert_every once overlapcheck yes all_in yes volumefraction_region 0.3 region bc
+ insert_every once overlapcheck yes all_in yes volumefraction_region 0.25 region bc
#fix ins nve_group pour/legacy 500 1 1 vol 0.7 1000 diam 0.008 0.008 dens 8000 8000 vel 0. 0. 0. 0. -0.3 region bc
diff --git a/examples/LIGGGHTS/Tutorials_public/heatTransfer_1/log.liggghts b/examples/LIGGGHTS/Tutorials_public/heatTransfer_1/log.liggghts
deleted file mode 100644
index 91b8d35..0000000
--- a/examples/LIGGGHTS/Tutorials_public/heatTransfer_1/log.liggghts
+++ /dev/null
@@ -1,202 +0,0 @@
-LIGGGHTS (Version LIGGGHTS-MASTER 2.3.5plusplus, compiled 2013-07-02-16:07:11 by ckloss based on LAMMPS 20 Apr 2012)
-#Heat transfer example
-
-atom_style granular
-atom_modify map array
-boundary m m m
-newton off
-
-communicate single vel yes
-
-units si
-
-region reg block -0.05 0.05 -0.05 0.05 0. 0.15 units box
-create_box 1 reg
-Created orthogonal box = (-0.05 -0.05 0) to (0.05 0.05 0.15)
- 1 by 1 by 1 MPI processor grid
-
-neighbor 0.002 bin
-neigh_modify delay 0
-
-#Material properties required for new pair styles
-
-fix m1 all property/global youngsModulus peratomtype 5.e6
-fix m2 all property/global poissonsRatio peratomtype 0.45
-fix m3 all property/global coefficientRestitution peratomtypepair 1 0.7
-fix m4 all property/global coefficientFriction peratomtypepair 1 0.05
-#fix m5 all property/global characteristicVelocity scalar 2.
-
-
-#New pair style
-pair_style gran/hertz/history #Hertzian without cohesion
-pair_coeff * *
-
-timestep 0.000025
-
-fix gravi all gravity 9.81 vector 0.0 0.0 -1.0
-
-fix zwalls1 all wall/gran/hertz/history primitive type 1 zplane 0.0
-fix zwalls2 all wall/gran/hertz/history primitive type 1 zplane 0.15
-fix cylwalls all wall/gran/hertz/history primitive type 1 zcylinder 0.05 0. 0.
-
-#heat transfer
-fix ftco all property/global thermalConductivity peratomtype 100.
-fix ftca all property/global thermalCapacity peratomtype 10.
-fix heattransfer all heat/gran initial_temperature 300.
-
-#region of insertion
-region bc cylinder z 0. 0. 0.045 0.00 0.15 units box
-
-#particle distributions and insertion
-fix pts1 all particletemplate/sphere 1 atom_type 1 density constant 8000 radius constant 0.004
-fix pdd1 all particledistribution/discrete 1. 1 pts1 1.0
-
-fix ins all insert/pack seed 100001 distributiontemplate pdd1 vel constant 0. 0. -0.3 insert_every once overlapcheck yes all_in yes volumefraction_region 0.3 region bc
-
-#fix ins nve_group pour/legacy 500 1 1 vol 0.7 1000 diam 0.008 0.008 dens 8000 8000 vel 0. 0. 0. 0. -0.3 region bc
-
-#apply nve integration to all particles
-fix integr all nve/sphere
-
-#output settings, include total thermal energy
-compute rke all erotate/sphere
-thermo_style custom step atoms ke c_rke f_heattransfer vol
-thermo 1000
-thermo_modify lost ignore norm no
-compute_modify thermo_temp dynamic yes
-
-#insert the first particles so that dump is not empty
-run 1
-INFO: Particle insertion ins: inserting every 0 steps
-Memory usage per processor = 10.246 Mbytes
-Step Atoms KinEng rke heattran Volume
- 0 0 -0 0 0 0.0015
-INFO: Particle insertion ins: inserted 838 particle templates (mass 1.797226) at step 1
- - a total of 838 particle templates (mass 1.797226) inserted so far.
- 1 838 0.08094128 0 5391.6767 0.0015
-Loop time of 0.0168691 on 1 procs for 1 steps with 838 atoms
-
-Pair time (%) = 2.19345e-05 (0.130028)
-Neigh time (%) = 0.000792027 (4.69514)
-Comm time (%) = 8.10623e-06 (0.0480538)
-Outpt time (%) = 2.5034e-05 (0.148402)
-Other time (%) = 0.016022 (94.9784)
-
-Nlocal: 838 ave 838 max 838 min
-Histogram: 1 0 0 0 0 0 0 0 0 0
-Nghost: 0 ave 0 max 0 min
-Histogram: 1 0 0 0 0 0 0 0 0 0
-Neighs: 1191 ave 1191 max 1191 min
-Histogram: 1 0 0 0 0 0 0 0 0 0
-
-Total # of neighbors = 1191
-Ave neighs/atom = 1.42124
-Neighbor list builds = 1
-Dangerous builds = 0
-dump dmp all custom 800 post/dump*.heatGran id type type x y z ix iy iz vx vy vz fx fy fz omegax omegay omegaz radius f_Temp[0] f_heatFlux[0]
-unfix ins
-
-#let the particles settle
-run 20000 upto
-Memory usage per processor = 11.0885 Mbytes
-Step Atoms KinEng rke heattran Volume
- 1 838 0.08094128 0 5391.6767 0.0015
- 1000 838 0.23330023 0.00028085647 5391.6767 0.0015
- 2000 838 0.40500607 0.00093155487 5391.6767 0.0015
- 3000 838 0.48121995 0.0020863291 5391.6767 0.0015
- 4000 838 0.30943228 0.0038350581 5391.6767 0.0015
- 5000 838 0.02882825 0.003107303 5391.6767 0.0015
- 6000 838 0.0065437378 0.0017489218 5391.6767 0.0015
- 7000 838 0.0053079981 0.0011778947 5391.6767 0.0015
- 8000 838 0.0039618257 0.00070482328 5391.6767 0.0015
- 9000 838 0.0021466144 0.00044086216 5391.6767 0.0015
- 10000 838 0.00083850394 0.0002768789 5391.6767 0.0015
- 11000 838 0.00039670374 0.00016310802 5391.6767 0.0015
- 12000 838 0.0001877567 9.8977974e-05 5391.6767 0.0015
- 13000 838 0.00010847408 5.4983425e-05 5391.6767 0.0015
- 14000 838 4.1837382e-05 2.5205873e-05 5391.6767 0.0015
- 15000 838 1.5219879e-05 1.2108599e-05 5391.6767 0.0015
- 16000 838 4.7182824e-06 5.1401193e-06 5391.6767 0.0015
- 17000 838 5.4179212e-06 3.1806858e-06 5391.6767 0.0015
- 18000 838 1.292492e-05 5.445621e-06 5391.6767 0.0015
- 19000 838 5.2689388e-07 1.473221e-06 5391.6767 0.0015
- 20000 838 3.2041393e-07 4.3607005e-07 5391.6767 0.0015
-Loop time of 10.8166 on 1 procs for 19999 steps with 838 atoms
-
-Pair time (%) = 6.72001 (62.1269)
-Neigh time (%) = 0.157796 (1.45884)
-Comm time (%) = 0.0115404 (0.106692)
-Outpt time (%) = 0.137965 (1.2755)
-Other time (%) = 3.78928 (35.0321)
-
-Nlocal: 838 ave 838 max 838 min
-Histogram: 1 0 0 0 0 0 0 0 0 0
-Nghost: 0 ave 0 max 0 min
-Histogram: 1 0 0 0 0 0 0 0 0 0
-Neighs: 3840 ave 3840 max 3840 min
-Histogram: 1 0 0 0 0 0 0 0 0 0
-
-Total # of neighbors = 3840
-Ave neighs/atom = 4.58234
-Neighbor list builds = 175
-Dangerous builds = 0
-
-#set particle temperature for half the bed
-region halfbed block 0 INF INF INF INF INF units box
-set region halfbed property/atom Temp 800.
- 424 settings made for property/atom
-
-#run to see heat transfer
-run 50000 upto
-Memory usage per processor = 11.0885 Mbytes
-Step Atoms KinEng rke heattran Volume
- 20000 838 3.2041393e-07 4.3607005e-07 9938.3572 0.0015
- 21000 838 1.6334093e-07 1.8420005e-07 9938.3572 0.0015
- 22000 838 7.1347443e-08 7.9822898e-08 9938.3572 0.0015
- 23000 838 8.9232552e-09 3.8361569e-08 9938.3572 0.0015
- 24000 838 6.2585103e-09 3.6690958e-08 9938.3572 0.0015
- 25000 838 4.8203317e-09 3.6041431e-08 9938.3572 0.0015
- 26000 838 3.9482728e-09 3.5895974e-08 9938.3572 0.0015
- 27000 838 3.5593562e-09 3.5712658e-08 9938.3572 0.0015
- 28000 838 3.4056738e-09 3.5624928e-08 9938.3572 0.0015
- 29000 838 3.2918842e-09 3.5639842e-08 9938.3572 0.0015
- 30000 838 2.4644172e-09 3.2981097e-08 9938.3572 0.0015
- 31000 838 2.2406862e-09 3.2940858e-08 9938.3572 0.0015
- 32000 838 6.3246241e-10 3.0695035e-08 9938.3572 0.0015
- 33000 838 1.7037387e-09 3.0063433e-08 9938.3572 0.0015
- 34000 838 1.2117871e-09 2.7510783e-08 9938.3572 0.0015
- 35000 838 1.1242344e-09 2.751392e-08 9938.3572 0.0015
- 36000 838 1.084762e-09 2.7512936e-08 9938.3572 0.0015
- 37000 838 1.066695e-09 2.7510646e-08 9938.3572 0.0015
- 38000 838 7.8090138e-10 2.6928417e-08 9938.3572 0.0015
- 39000 838 8.4800071e-10 2.6595997e-08 9938.3572 0.0015
- 40000 838 8.5699603e-10 2.6595621e-08 9938.3572 0.0015
- 41000 838 8.6207934e-10 2.6595459e-08 9938.3572 0.0015
- 42000 838 8.5929449e-10 2.6595407e-08 9938.3572 0.0015
- 43000 838 8.5162952e-10 2.6595109e-08 9938.3572 0.0015
- 44000 838 8.4408513e-10 2.6594904e-08 9938.3572 0.0015
- 45000 838 8.3830893e-10 2.6594733e-08 9938.3572 0.0015
- 46000 838 8.3321612e-10 2.6594399e-08 9938.3572 0.0015
- 47000 838 8.2126987e-10 2.6588734e-08 9938.3572 0.0015
- 48000 838 5.861104e-10 2.5010593e-08 9938.3572 0.0015
- 49000 838 5.6959467e-10 2.4964372e-08 9938.3572 0.0015
- 50000 838 5.6903337e-10 2.4964193e-08 9938.3572 0.0015
-Loop time of 19.1494 on 1 procs for 30000 steps with 838 atoms
-
-Pair time (%) = 11.9029 (62.1581)
-Neigh time (%) = 0 (0)
-Comm time (%) = 0.01162 (0.060681)
-Outpt time (%) = 0.605505 (3.16201)
-Other time (%) = 6.62937 (34.6192)
-
-Nlocal: 838 ave 838 max 838 min
-Histogram: 1 0 0 0 0 0 0 0 0 0
-Nghost: 0 ave 0 max 0 min
-Histogram: 1 0 0 0 0 0 0 0 0 0
-Neighs: 3841 ave 3841 max 3841 min
-Histogram: 1 0 0 0 0 0 0 0 0 0
-
-Total # of neighbors = 3841
-Ave neighs/atom = 4.58353
-Neighbor list builds = 0
-Dangerous builds = 0
diff --git a/examples/LIGGGHTS/Tutorials_public/heatTransfer_2/in.heatGran b/examples/LIGGGHTS/Tutorials_public/heatTransfer_2/in.heatGran
index 20b453d..216a11b 100644
--- a/examples/LIGGGHTS/Tutorials_public/heatTransfer_2/in.heatGran
+++ b/examples/LIGGGHTS/Tutorials_public/heatTransfer_2/in.heatGran
@@ -49,7 +49,7 @@ fix pts1 all particletemplate/sphere 1 atom_type 1 density constant 8000 radius
fix pdd1 all particledistribution/discrete 1. 1 pts1 1.0
fix ins all insert/pack seed 100001 distributiontemplate pdd1 vel constant 0. 0. -0.3 &
- insert_every once overlapcheck yes all_in yes volumefraction_region 0.3 region bc
+ insert_every once overlapcheck yes all_in yes volumefraction_region 0.25 region bc
#fix ins nve_group pour/legacy 500 1 1 vol 0.7 1000 diam 0.008 0.008 dens 8000 8000 vel 0. 0. 0. 0. -0.3 region bc
diff --git a/src/atom_vec_ellipsoid.cpp b/src/atom_vec_ellipsoid.cpp
index fc0396e..9dae26b 100644
--- a/src/atom_vec_ellipsoid.cpp
+++ b/src/atom_vec_ellipsoid.cpp
@@ -656,6 +656,7 @@ int AtomVecEllipsoid::pack_border_vel(int n, int *list, double *buf,
if (ellipsoid[j] < 0) buf[m++] = 0;
else {
buf[m++] = 1;
+ shape = bonus[ellipsoid[j]].shape;
quat = bonus[ellipsoid[j]].quat;
buf[m++] = shape[0];
buf[m++] = shape[1];
diff --git a/src/atom_vec_sphere.cpp b/src/atom_vec_sphere.cpp
index a12cfae..811fbb2 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 a908835..a73085c 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/domain_wedge_dummy.h b/src/atom_vec_sphere_w.cpp
similarity index 62%
copy from src/domain_wedge_dummy.h
copy to src/atom_vec_sphere_w.cpp
index 8a6bee9..adffed4 100644
--- a/src/domain_wedge_dummy.h
+++ b/src/atom_vec_sphere_w.cpp
@@ -19,29 +19,28 @@
See the README file in the top-level directory.
------------------------------------------------------------------------- */
-/* ----------------------------------------------------------------------
- Contributing authors:
- Stefan Amberger (JKU Linz)
- Christoph Kloss (JKU Linz, DCS Computing GmbH, Linz)
-------------------------------------------------------------------------- */
+#include "atom_vec_sphere.h"
+#include "domain_wedge.h"
-#ifndef DOMAIN_WEDGE_H
-#define DOMAIN_WEDGE_H
+#ifndef DOMAIN_WEDGE_REAL_H
+#define DOMAIN_WEDGE_REAL_H
-#include "domain.h"
+using namespace LAMMPS_NS;
-namespace LAMMPS_NS {
+/* ---------------------------------------------------------------------- */
-class DomainWedge : public Domain
+int AtomVecSphere::pack_border_vel_wedge(int n, int *list, double *buf,
+ int pbc_flag, int *pbc)
{
+ return 0;
+}
- public:
-
- DomainWedge(class LAMMPS *lmp) : Domain(lmp) {};
- void set_domain(class RegWedge *rw) {}
-
-};
+/* ---------------------------------------------------------------------- */
+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 619c950..94cae84 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 586cc5d..88759e5 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 0000000..45afbd8
--- /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 at 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 at 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 f0ba88c..20e9a09 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 815a665..3e41f30 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 3efd698..6592730 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 5dbfeb3..88406c0 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 5709c42..869d56e 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 3c815ed..ccd91c9 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/compute_reduce.cpp b/src/compute_reduce.cpp
index 645dbac..1736cba 100644
--- a/src/compute_reduce.cpp
+++ b/src/compute_reduce.cpp
@@ -30,7 +30,7 @@
using namespace LAMMPS_NS;
enum{SUM,MINN,MAXX,AVE};
-enum{X,V,F,COMPUTE,FIX,VARIABLE};
+enum{X,V,F,RHO,P,COMPUTE,FIX,VARIABLE};
enum{PERATOM,LOCAL};
#define INVOKED_VECTOR 2
@@ -111,6 +111,12 @@ ComputeReduce::ComputeReduce(LAMMPS *lmp, int narg, char **arg) :
} else if (strcmp(arg[iarg],"fz") == 0) {
which[nvalues] = F;
argindex[nvalues++] = 2;
+ } else if (strcmp(arg[iarg],"rho") == 0) {
+ which[nvalues] = RHO;
+ argindex[nvalues++] = 0;
+ } else if (strcmp(arg[iarg],"p") == 0) {
+ which[nvalues] = P;
+ argindex[nvalues++] = 0;
} else if (strncmp(arg[iarg],"c_",2) == 0 ||
strncmp(arg[iarg],"f_",2) == 0 ||
@@ -471,6 +477,18 @@ double ComputeReduce::compute_one(int m, int flag)
for (i = 0; i < nlocal; i++)
if (mask[i] & groupbit) combine(one,f[i][aidx],i);
} else one = f[flag][aidx];
+ } else if (which[m] == RHO) {
+ double *rho = atom->rho;
+ if (flag < 0) {
+ for (i = 0; i < nlocal; i++)
+ if (mask[i] & groupbit) combine(one,rho[i],i);
+ } else one = rho[flag];
+ } else if (which[m] == P) {
+ double *p = atom->p;
+ if (flag < 0) {
+ for (i = 0; i < nlocal; i++)
+ if (mask[i] & groupbit) combine(one,p[i],i);
+ } else one = p[flag];
// invoke compute if not previously invoked
diff --git a/src/container_base.cpp b/src/container_base.cpp
index ca514fb..4ce62a0 100644
--- a/src/container_base.cpp
+++ b/src/container_base.cpp
@@ -45,11 +45,23 @@ using namespace LAMMPS_NS;
{
}
+ ContainerBase::ContainerBase(char *_id)
+ : id_(0),
+ communicationType_(COMM_TYPE_MANUAL),
+ refFrame_(REF_FRAME_UNDEFINED),
+ restartType_(RESTART_TYPE_UNDEFINED),
+ scalePower_(-1)
+ {
+ id_ = new char[strlen(_id)+1];
+ strcpy(id_,_id);
+ }
+
ContainerBase::ContainerBase(char *_id, char* _comm, char* _ref, char *_restart,int _scalePower)
: id_(0),
communicationType_(COMM_TYPE_MANUAL),
+ refFrame_(REF_FRAME_UNDEFINED),
restartType_(RESTART_TYPE_UNDEFINED),
- refFrame_(REF_FRAME_UNDEFINED)
+ scalePower_(-1)
{
setProperties(_id, _comm, _ref,_restart,_scalePower);
}
@@ -66,7 +78,7 @@ using namespace LAMMPS_NS;
ContainerBase::~ContainerBase()
{
- delete []id_;
+ if(id_) delete []id_;
}
/* ----------------------------------------------------------------------
diff --git a/src/container_base.h b/src/container_base.h
index 3bddf24..bed638f 100644
--- a/src/container_base.h
+++ b/src/container_base.h
@@ -45,6 +45,8 @@ namespace LAMMPS_NS
{
public:
+ ContainerBase(char *_id);
+
virtual ~ContainerBase();
void setProperties(char *_id, char* _comm, char* _ref, char *_restart,int _scalePower = 1);
@@ -105,7 +107,6 @@ namespace LAMMPS_NS
protected:
- ContainerBase();
ContainerBase(char *_id, char* _comm, char* _ref, char *_restart,int _scalePower);
ContainerBase(ContainerBase const &orig);
@@ -124,6 +125,10 @@ namespace LAMMPS_NS
int refFrame_;
int scalePower_;
int restartType_;
+
+ private:
+
+ ContainerBase();
};
// *************************************
diff --git a/src/container_base_I.h b/src/container_base_I.h
index d194d41..27f2c8f 100644
--- a/src/container_base_I.h
+++ b/src/container_base_I.h
@@ -80,33 +80,33 @@
{
// return true for manual communication, such as for node_, node_orig_
// etc in MultiNodeMeshParallel
- if(communicationType_ == COMM_TYPE_MANUAL)
+ if(COMM_TYPE_MANUAL == communicationType_)
return true;
- if(operation == OPERATION_RESTART)
+ if(OPERATION_RESTART == operation)
{
if(restartType_ == RESTART_TYPE_YES)
return true;
return false;
}
- if(operation == OPERATION_COMM_BORDERS ||
- operation == OPERATION_COMM_EXCHANGE )
+ if(OPERATION_COMM_BORDERS == operation ||
+ OPERATION_COMM_EXCHANGE == operation )
return true;
- if(communicationType_ == COMM_TYPE_NONE)
+ if(COMM_TYPE_NONE == communicationType_)
return false;
- if(operation == OPERATION_COMM_REVERSE &&
- communicationType_ == COMM_TYPE_REVERSE)
+ if(OPERATION_COMM_REVERSE == operation &&
+ COMM_TYPE_REVERSE == communicationType_)
return true;
- if(operation == OPERATION_COMM_FORWARD &&
- communicationType_ == COMM_TYPE_FORWARD)
+ if(OPERATION_COMM_FORWARD == operation &&
+ COMM_TYPE_FORWARD == communicationType_)
return true;
- if(operation == OPERATION_COMM_FORWARD &&
- communicationType_ == COMM_TYPE_FORWARD_FROM_FRAME)
+ if(OPERATION_COMM_FORWARD == operation &&
+ COMM_TYPE_FORWARD_FROM_FRAME == communicationType_)
{
if(scale && !isScaleInvariant())
return true;
diff --git a/src/debug_liggghts.h b/src/debug_liggghts.h
index 4f6430b..3d01933 100644
--- a/src/debug_liggghts.h
+++ b/src/debug_liggghts.h
@@ -28,6 +28,7 @@
#include "style_fix.h"
#include "vector_liggghts.h"
#include "container.h"
+#include "atom.h"
namespace LAMMPS_NS {
diff --git a/src/domain.cpp b/src/domain.cpp
index fd72abc..d73e79f 100644
--- a/src/domain.cpp
+++ b/src/domain.cpp
@@ -44,6 +44,7 @@
#include "math_const.h"
#include "memory.h"
#include "error.h"
+#include "neighbor.h"
using namespace LAMMPS_NS;
using namespace MathConst;
@@ -93,6 +94,7 @@ Domain::Domain(LAMMPS *lmp) : Pointers(lmp)
nregion = maxregion = 0;
regions = NULL;
+ is_wedge = false;
}
/* ---------------------------------------------------------------------- */
@@ -1415,40 +1417,3 @@ void Domain::box_corners()
lamda2x(corners[7],corners[7]);
}
-/* ----------------------------------------------------------------------
- domain checks - 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];
-
- if(i < nlocal) return 0;
-
- else
- {
- for(idim = 0; idim < 3; idim++)
- if ((x[idim] < boxlo[idim] || x[idim] > boxhi[idim]) && periodicity[idim])
- return 1;
- }
- return 0;
-}
-
-int Domain::is_periodic_ghost_of_owned(int i)
-{
- int idim;
- int nlocal = atom->nlocal;
- double *x = atom->x[i];
-
- if(i < nlocal) return 0;
-
- else
- {
- for(idim = 0; idim < 3; idim++)
- if ((x[idim] < boxlo[idim] || x[idim] > boxhi[idim]) && periodicity[idim] && 1 == comm->procgrid[idim])
- return 1;
- }
- return 0;
-}
diff --git a/src/domain.h b/src/domain.h
index 8c5abea..fb4f37d 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);
- int is_periodic_ghost_of_owned(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 0000000..54bc314
--- /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 at 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 at 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 8a6bee9..6026d61 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/dump_mesh_vtk.cpp b/src/dump_mesh_vtk.cpp
index 30a581e..3dbee03 100644
--- a/src/dump_mesh_vtk.cpp
+++ b/src/dump_mesh_vtk.cpp
@@ -667,8 +667,8 @@ void DumpMeshVTK::write_data_ascii_point(int n, double *mybuf)
m = 0;
buf_pos = 0;
- ScalarContainer<double> points;
- ScalarContainer<int> tri_points;
+ ScalarContainer<double> points("DumpMeshVTK::points");
+ ScalarContainer<int> tri_points("DumpMeshVTK::tri_points");
bool add;
for (int i = 0; i < n; i++)
{
@@ -703,7 +703,7 @@ void DumpMeshVTK::write_data_ascii_point(int n, double *mybuf)
class ScalarContainer<int> **points_neightri;
points_neightri = new ScalarContainer<int>*[(int)points.size()/3];
for (int i=0; i < points.size()/3; i++)
- points_neightri[i] = new ScalarContainer<int>;
+ points_neightri[i] = new ScalarContainer<int>("DumpMeshVTK::points_neightri");
for (int i=0; i < 3*n; i+=3)
{
for (int j=0; j<3;j++)
diff --git a/src/fix.h b/src/fix.h
index 174b9ec..4d3b27a 100644
--- a/src/fix.h
+++ b/src/fix.h
@@ -98,9 +98,12 @@ class Fix : protected Pointers {
virtual int setmask() = 0;
+ 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_ave_spatial.cpp b/src/fix_ave_spatial.cpp
index d69e7cc..734cfa9 100644
--- a/src/fix_ave_spatial.cpp
+++ b/src/fix_ave_spatial.cpp
@@ -337,8 +337,8 @@ FixAveSpatial::FixAveSpatial(LAMMPS *lmp, int narg, char **arg) :
if (fp2 && me == 0) {
fprintf(fp2,"# Spatial-averaged data for fix %s and group %s\n",
id,arg[1]);
- fprintf(fp2,"# Mean and standard deviation for all bins including between L=%i and H=%i particles\n",lowerLimit,upperLimit);
- fprintf(fp2,"# Timestep Natoms Nbins MaxAtomsPerBin NbinsEmpty Nbins<L Nbins>H Nsamples NatomsPerSample ");
+ fprintf(fp2,"# Mean and standard deviation for all bins including between N1=%i and N2=%i particles\n",lowerLimit,upperLimit);
+ fprintf(fp2,"# Timestep Natoms Nbins MaxAtomsPerBin NbinsEmpty Nbins<N1 Nbins>N2 Nsamples NatomsPerSample ");
for (int i = 0; i < nvalues; i++) {
fprintf(fp2,"{%s: trueMean samplesMean Std} ",arg[6+3*ndim+i]);
@@ -880,7 +880,7 @@ void FixAveSpatial::end_of_step()
if (i == 0) count_sum += count_total[m];
}
- true_mean[i] /= count_sum;
+ if (count_sum != 0) true_mean[i] /= count_sum;
samples_mean[i] = 0;
std[i] = 0;
diff --git a/src/fix_cfd_coupling_force.cpp b/src/fix_cfd_coupling_force.cpp
index eefec91..6b5326f 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 49e41c3..1384f4f 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;
+
}
/* ---------------------------------------------------------------------- */
@@ -55,6 +79,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 e4918fc..7beddfa 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 5fa7f8a..301f2f9 100644
--- a/src/fix_contact_history.cpp
+++ b/src/fix_contact_history.cpp
@@ -71,7 +71,9 @@ FixContactHistory::FixContactHistory(LAMMPS *lmp, int narg, char **arg) :
//read dnum
dnum = atoi(arg[iarg++]);
-
+
+ index_decide_noncontacting = -1;
+
if(dnum < 0)
error->all(FLERR,"dnum must be >=0 in fix contacthistory");
if(dnum > 10)
@@ -191,6 +193,9 @@ void FixContactHistory::init()
if (atom->tag_enable == 0)
error->fix_error(FLERR,this,"using contact history requires atoms have IDs");
+ if (-1 == index_decide_noncontacting && 1. < neighbor->contactHistoryDistanceFactor)
+ error->fix_error(FLERR,this,"have to call set_index_decide_noncontacting() function");
+
if(is_pair)
{
if(!force->pair_match("gran", 0))
@@ -209,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)
@@ -286,7 +291,7 @@ void FixContactHistory::pre_exchange_pair()
int *tag = atom->tag;
double **x = atom->x;
- double *radius = atom->radius;
+ double *radius = atom->radius;
NeighList *list = pair_gran->list;
inum = list->inum;
@@ -295,7 +300,7 @@ void FixContactHistory::pre_exchange_pair()
firstneigh = list->firstneigh;
firsttouch = list->listgranhistory->firstneigh;
firsthist = list->listgranhistory->firstdouble;
- contactHistDistanceFactor = neighbor->contactHistoryDistanceFactor;
+ contactHistDistanceFactor = neighbor->contactHistoryDistanceFactor;
if(contactHistDistanceFactor> 1.0) considerNonContactingParticles = true;
for (ii = 0; ii < inum; ii++) {
@@ -312,28 +317,20 @@ void FixContactHistory::pre_exchange_pair()
for (jj = 0; jj < jnum; jj++) {
j = jlist[jj];
j &= NEIGHMASK;
-#if 0
+
//Check if considerNonContactingParticles are within range
- haveNonContactingParticlesInRange = false;
- if(considerNonContactingParticles && (j<nlocal))
- {
- radSum = radius[j] + rPartner;
- delx = x[j][0] - xPartner[0];
- dely = x[j][1] - xPartner[1];
- delz = x[j][2] - xPartner[2];
- rsq = delx*delx + dely*dely + delz*delz;
- if( rsq<(contactHistDistanceFactor*radSum*contactHistDistanceFactor*radSum) ) //check if particles are close enough to keep contact history
- haveNonContactingParticlesInRange = true;
- }
+ haveNonContactingParticlesInRange = false;
+
+ if(considerNonContactingParticles)
+ {
+ hist = &allhist[dnum*jj];
+ if( hist[index_decide_noncontacting]>0.0 ) //check if particles are close enough to keep contact history
+ haveNonContactingParticlesInRange = true;
+ }
-// fprintf(screen, "***FixContactHistory::pre_exchange_pair - haveNonContactingParticlesInRange %d, considerNonContactingParticles %d, contactHistDistanceFactor: %f \n.",
-// haveNonContactingParticlesInRange, considerNonContactingParticles, contactHistDistanceFactor);
-// Check if we need to consider non-contacting particles
if ( (touch[jj] ) || haveNonContactingParticlesInRange)
-#endif
- if ( (touch[jj] ) || considerNonContactingParticles)
- {
- hist = &allhist[dnum*jj];
+ {
+ hist = &allhist[dnum*jj];
if (npartner[i] >= maxtouch) grow_arrays_maxtouch(atom->nmax);
m = npartner[i];
diff --git a/src/fix_contact_history.h b/src/fix_contact_history.h
index a3140ed..b8bf57c 100644
--- a/src/fix_contact_history.h
+++ b/src/fix_contact_history.h
@@ -81,6 +81,9 @@ class FixContactHistory : public Fix {
int n_contacts();
int n_contacts(int contact_groupbit);
+ void set_index_decide_noncontacting(int index)
+ { index_decide_noncontacting = index; };
+
int get_dnum()
{ return dnum; }
@@ -118,6 +121,9 @@ class FixContactHistory : public Fix {
int dnum;
int *newtonflag;
char **history_id;
+
+ //
+ int index_decide_noncontacting;
};
// *************************************
diff --git a/src/fix_heat_gran.cpp b/src/fix_heat_gran.cpp
index f43fd78..ceb995e 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 5d7475b..9dce5dd 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. "
@@ -450,9 +451,6 @@ void FixInsert::init()
if(!fix_multisphere) multisphere = NULL;
else multisphere = &fix_multisphere->data();
- if(fix_multisphere && fix_multisphere->igroup != igroup)
- error->fix_error(FLERR,this,"Fix insert command and fix multisphere command are not compatible, must be same group");
-
// in case of new fix insert in a restarted simulation, have to add current time-step
if(next_reneighbor > 0 && next_reneighbor < ntimestep)
error->fix_error(FLERR,this,"'start' step can not be before current step");
diff --git a/src/fix_insert_pack.cpp b/src/fix_insert_pack.cpp
index 897a3a9..de63343 100644
--- a/src/fix_insert_pack.cpp
+++ b/src/fix_insert_pack.cpp
@@ -150,7 +150,7 @@ void FixInsertPack::calc_insertion_properties()
ins_region->reset_random(seed + SEED_OFFSET);
calc_region_volume_local();
- if(region_volume <= 0. || region_volume_local < 0. || region_volume_local > region_volume)
+ if(region_volume <= 0. || region_volume_local < 0. || (region_volume_local - region_volume)/region_volume > 1e-3 )
error->one(FLERR,"Fix insert: Region volume calculation with MC failed");
if(ins_region->dynamic_check())
@@ -187,6 +187,7 @@ void FixInsertPack::calc_region_volume_local()
{
ins_region->volume_mc(ntry_mc,all_in_flag==0?false:true,fix_distribution->max_r_bound(),
region_volume,region_volume_local);
+
}
/* ----------------------------------------------------------------------
diff --git a/src/fix_massflow_mesh.cpp b/src/fix_massflow_mesh.cpp
index 982eb26..83c0443 100644
--- a/src/fix_massflow_mesh.cpp
+++ b/src/fix_massflow_mesh.cpp
@@ -23,6 +23,7 @@
#include "stdlib.h"
#include "string.h"
#include "atom.h"
+#include "atom_vec.h"
#include "comm.h"
#include "modify.h"
#include "force.h"
@@ -55,10 +56,14 @@ FixMassflowMesh::FixMassflowMesh(LAMMPS *lmp, int narg, char **arg) :
mass_last_(0.),
t_count_(0.),
delta_t_(0.),
+ fp_(0),
reset_t_count_(true),
havePointAtOutlet_(false),
insideOut_(false),
- screenflag_(false)
+ screenflag_(false),
+ delete_atoms_(false),
+ mass_deleted_(0.),
+ nparticles_deleted_(0)
{
vectorZeroize3D(nvec_);
vectorZeroize3D(pref_);
@@ -67,10 +72,6 @@ FixMassflowMesh::FixMassflowMesh(LAMMPS *lmp, int narg, char **arg) :
// parse args for this class
- if(1 < comm->nprocs && 0 == comm->me)
- fprintf(screen,"**FixMassflowMesh: > 1 process - "
- " will write to multiple files\n");
-
int iarg = 3;
bool hasargs = true;
@@ -153,14 +154,28 @@ FixMassflowMesh::FixMassflowMesh(LAMMPS *lmp, int narg, char **arg) :
else error->all(FLERR,"Illegal fix print command");
iarg += 2;
hasargs = true;
+ } else if (strcmp(arg[iarg],"delete_atoms") == 0) {
+ if(narg < iarg+2)
+ error->fix_error(FLERR,this,"Illegal keyword entry");
+ if (strcmp(arg[iarg+1],"yes") == 0) delete_atoms_ = true;
+ else if (strcmp(arg[iarg+1],"no") == 0) delete_atoms_ = false;
+ else error->all(FLERR,"Illegal delete command");
+ iarg += 2;
+ hasargs = true;
} else
error->fix_error(FLERR,this,"unknown keyword");
}
+ if(fp_ && 1 < comm->nprocs && 0 == comm->me)
+ fprintf(screen,"**FixMassflowMesh: > 1 process - "
+ " will write to multiple files\n");
+
// error checks on necessary args
+ if(!once_ && delete_atoms_)
+ error->fix_error(FLERR,this,"using 'delete_atoms yes' requires 'count once'");
if( !once_ && havePointAtOutlet_)
- error->fix_error(FLERR,this,"setting 'pointAtOutlet' requires 'count once'");
+ error->fix_error(FLERR,this,"setting 'point_at_outlet' requires 'count once'");
if( vectorMag3D(sidevec_)==0. && !havePointAtOutlet_)
error->fix_error(FLERR,this,"expecting keyword 'vec_side'");
if (!fix_mesh_)
@@ -240,6 +255,9 @@ void FixMassflowMesh::init()
if(modify->n_fixes_style("multisphere"))
error->fix_error(FLERR,this,"does not support multi-sphere");
+
+ if(delete_atoms_ && 0 == atom->map_style)
+ error->fix_error(FLERR,this,"requires an 'atom_modify map' command to allocate an atom map");
}
/* ---------------------------------------------------------------------- */
@@ -248,8 +266,8 @@ void FixMassflowMesh::setup(int vflag)
{
// check if face planar
- if(!fix_mesh_->triMesh()->isPlanar())
- error->fix_error(FLERR,this,"command requires a planar face mass flow measurement");
+ if(!fix_mesh_->triMesh()->isPlanar() && !havePointAtOutlet_)
+ error->fix_error(FLERR,this,"requires a planar face mass flow measurement or using 'point_at_outlet'");
}
/* ---------------------------------------------------------------------- */
@@ -258,6 +276,7 @@ int FixMassflowMesh::setmask()
{
int mask = 0;
mask |= POST_INTEGRATE;
+ if(delete_atoms_) mask |= PRE_EXCHANGE;
return mask;
}
@@ -328,9 +347,7 @@ void FixMassflowMesh::post_integrate()
dot = vectorDot3D(delta,nvec_);
}
else //particle is not overlapping with mesh, so continue
- {
continue;
- }
if(insideOut_) dot =-dot;
}
@@ -354,6 +371,9 @@ void FixMassflowMesh::post_integrate()
{
mass_this += rmass[iPart];
nparticles_this ++;
+ if(delete_atoms_)
+ atom_tags_delete_.push_back(atom->tag[iPart]);
+
if (screenflag_ && screen)
fprintf(screen," %d %4.4g %4.4g %4.4g %4.4g %4.4g %4.4g %4.4g \n ",
tag[iPart],2.*radius[iPart]/force->cg(),
@@ -388,17 +408,81 @@ void FixMassflowMesh::post_integrate()
}
/* ----------------------------------------------------------------------
+ perform particle deletion of marked particles
+ done before exchange, borders, reneighbor
+ so that ghost atoms and neighbor lists will be correct
+------------------------------------------------------------------------- */
+
+void FixMassflowMesh::pre_exchange()
+{
+ int nlocal = atom->nlocal;
+ int iPart;
+ double *counter = fix_counter_->vector_atom;
+
+ if (update->ntimestep != next_reneighbor) return;
+
+ if (delete_atoms_)
+ {
+ double mass_deleted_this_ = 0.;
+ int nparticles_deleted_this_ = 0.;
+
+ // delete particles
+
+ while (atom_tags_delete_.size() > 0)
+ {
+ iPart = atom->map(atom_tags_delete_[0]);
+
+ mass_deleted_this_ += atom->rmass[iPart];
+ nparticles_deleted_this_++;
+
+ atom->avec->copy(atom->nlocal-1,iPart,1);
+ atom->nlocal--;
+
+ atom_tags_delete_.erase(atom_tags_delete_.begin());
+ }
+
+ atom_tags_delete_.clear();
+
+ MPI_Sum_Scalar(mass_deleted_this_,world);
+ MPI_Sum_Scalar(nparticles_deleted_this_,world);
+
+ mass_deleted_ += mass_deleted_this_;
+ nparticles_deleted_ += nparticles_deleted_this_;
+
+ if(nparticles_deleted_this_)
+ {
+ int i;
+ if (atom->molecular == 0) {
+ int *tag = atom->tag;
+ for (i = 0; i < atom->nlocal; i++) tag[i] = 0;
+ atom->tag_extend();
+ }
+
+ if (atom->tag_enable) {
+ if (atom->map_style) {
+ atom->nghost = 0;
+ atom->map_init();
+ atom->map_set();
+ }
+ }
+ }
+ }
+}
+
+/* ----------------------------------------------------------------------
pack entire state of Fix into one write
------------------------------------------------------------------------- */
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_;
list[n++] = static_cast<double>(nparticles_last_);
+ list[n++] = mass_deleted_;
+ list[n++] = static_cast<double>(nparticles_deleted_);
if (comm->me == 0) {
int size = n * sizeof(double);
@@ -420,7 +504,8 @@ void FixMassflowMesh::restart(char *buf)
t_count_ = list[n++];
mass_last_ = list[n++];
nparticles_last_ = static_cast<int>(list[n++]);
- nparticles_last_ = static_cast<int>(list[n++]);
+ mass_deleted_ = list[n++];
+ nparticles_deleted_ = static_cast<int>(list[n++]);
}
/* ----------------------------------------------------------------------
@@ -444,6 +529,10 @@ double FixMassflowMesh::compute_vector(int index)
return delta_t_ == 0. ? 0. : (mass_-mass_last_)/delta_t_;
if(index == 3)
return delta_t_ == 0. ? 0. : static_cast<double>(nparticles_-nparticles_last_)/delta_t_;
+ if(index == 4)
+ return mass_deleted_;
+ if(index == 5)
+ return static_cast<double>(nparticles_deleted_);
return 0.;
}
diff --git a/src/fix_massflow_mesh.h b/src/fix_massflow_mesh.h
index 711ff69..10971e6 100644
--- a/src/fix_massflow_mesh.h
+++ b/src/fix_massflow_mesh.h
@@ -29,6 +29,9 @@ FixStyle(massflow/mesh,FixMassflowMesh)
#define LMP_FIX_MASSFLOW_MESH_H
#include "fix.h"
+#include <vector>
+
+using namespace std;
namespace LAMMPS_NS {
@@ -47,6 +50,7 @@ class FixMassflowMesh : public Fix {
int setmask();
void post_integrate();
+ void pre_exchange();
void write_restart(FILE *fp);
void restart(char *buf);
@@ -84,6 +88,12 @@ class FixMassflowMesh : public Fix {
double t_count_, delta_t_;
bool reset_t_count_;
+ // in case particles counted should be deleted
+ bool delete_atoms_;
+ vector<int> atom_tags_delete_;
+ double mass_deleted_;
+ double nparticles_deleted_;
+
}; //end class
}
diff --git a/src/fix_mesh.cpp b/src/fix_mesh.cpp
index cb64b3f..05f87c1 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 e84b942..a75f53d 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)
@@ -119,13 +119,10 @@ FixMeshSurfaceStress::~FixMeshSurfaceStress()
{
}
-
/* ---------------------------------------------------------------------- */
-void FixMeshSurfaceStress::post_create()
+void FixMeshSurfaceStress::post_create_pre_restart()
{
- FixMeshSurface::post_create();
-
if(stress_flag_)
initStress();
@@ -135,6 +132,13 @@ void FixMeshSurfaceStress::post_create()
/* ---------------------------------------------------------------------- */
+void FixMeshSurfaceStress::post_create()
+{
+ FixMeshSurface::post_create();
+}
+
+/* ---------------------------------------------------------------------- */
+
void FixMeshSurfaceStress::init()
{
if(stress_flag_)
@@ -224,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
@@ -261,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;
@@ -283,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.h b/src/fix_mesh_surface_stress.h
index 5ee7df8..0e5e548 100644
--- a/src/fix_mesh_surface_stress.h
+++ b/src/fix_mesh_surface_stress.h
@@ -45,6 +45,7 @@ namespace LAMMPS_NS
FixMeshSurfaceStress(LAMMPS *lmp, int narg, char **arg);
virtual ~FixMeshSurfaceStress();
+ virtual void post_create_pre_restart();
virtual void post_create();
virtual void init();
diff --git a/src/fix_mesh_surface_stress_servo.cpp b/src/fix_mesh_surface_stress_servo.cpp
index 457ec08..648564c 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./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./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./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./set_point_;
+ set_point_inv_ = 1./fabs(set_point_);
sp_style_ = CONSTANT;
}
diff --git a/src/fix_move.cpp b/src/fix_move.cpp
index dc146c7..2d7641a 100644
--- a/src/fix_move.cpp
+++ b/src/fix_move.cpp
@@ -748,6 +748,7 @@ void FixMove::final_integrate()
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
if (xflag)
+ {
if (rmass) {
dtfm = dtf / rmass[i];
v[i][0] += dtfm * f[i][0];
@@ -755,8 +756,10 @@ void FixMove::final_integrate()
dtfm = dtf / mass[type[i]];
v[i][0] += dtfm * f[i][0];
}
+ }
if (yflag)
+ {
if (rmass) {
dtfm = dtf / rmass[i];
v[i][1] += dtfm * f[i][1];
@@ -764,8 +767,10 @@ void FixMove::final_integrate()
dtfm = dtf / mass[type[i]];
v[i][1] += dtfm * f[i][1];
}
+ }
if (zflag)
+ {
if (rmass) {
dtfm = dtf / rmass[i];
v[i][2] += dtfm * f[i][2];
@@ -773,6 +778,7 @@ void FixMove::final_integrate()
dtfm = dtf / mass[type[i]];
v[i][2] += dtfm * f[i][2];
}
+ }
}
}
}
diff --git a/src/fix_neighlist_mesh.cpp b/src/fix_neighlist_mesh.cpp
index 87f6ff2..be4735b 100644
--- a/src/fix_neighlist_mesh.cpp
+++ b/src/fix_neighlist_mesh.cpp
@@ -43,6 +43,8 @@ using namespace FixConst;
FixNeighlistMesh::FixNeighlistMesh(LAMMPS *lmp, int narg, char **arg)
: Fix(lmp,narg,arg),
+ contactList("contactList"),
+ numContacts("numContacts"),
buildNeighList(false),
movingMesh(false)
{
diff --git a/src/fix_property_atom_tracer.cpp b/src/fix_property_atom_tracer.cpp
index 9a5c1e5..7ac92f7 100644
--- a/src/fix_property_atom_tracer.cpp
+++ b/src/fix_property_atom_tracer.cpp
@@ -95,8 +95,8 @@ FixPropertyAtomTracer::FixPropertyAtomTracer(LAMMPS *lmp, int narg, char **arg,b
error->fix_error(FLERR,this,"not enough arguments for 'mark_step'");
iarg_++;
step_ = atoi(arg[iarg_++]);
- if(step_ < 0)
- error->fix_error(FLERR,this,"mark_step > 0 required");
+ if(step_ < 0 || step_ < update->ntimestep)
+ error->fix_error(FLERR,this,"mark_step > 0 required, mark_step must not be before current time-step");
hasargs = true;
} else if(strcmp(arg[iarg_],"marker_style") == 0) {
if(narg < iarg_+2)
@@ -106,6 +106,8 @@ FixPropertyAtomTracer::FixPropertyAtomTracer(LAMMPS *lmp, int narg, char **arg,b
marker_style_ = MARKER_HEAVISIDE;
else if(strcmp(arg[iarg_],"dirac") == 0)
marker_style_ = MARKER_DIRAC;
+ else if(strcmp(arg[iarg_],"none") == 0)
+ marker_style_ = MARKER_NONE;
else
error->fix_error(FLERR,this,"expecting 'heaviside' or 'dirac' after keyword 'marker_style'");
iarg_++;
@@ -176,7 +178,7 @@ void FixPropertyAtomTracer::end_of_step()
{
int ts = update->ntimestep;
- if(ts < step_ || (marker_style_ == MARKER_DIRAC && !first_mark_))
+ if(ts < step_ || marker_style_ == MARKER_NONE || (marker_style_ == MARKER_DIRAC && !first_mark_))
return;
int nlocal = atom->nlocal;
diff --git a/src/fix_property_atom_tracer.h b/src/fix_property_atom_tracer.h
index e0f0dee..f806b06 100644
--- a/src/fix_property_atom_tracer.h
+++ b/src/fix_property_atom_tracer.h
@@ -35,7 +35,8 @@ namespace LAMMPS_NS {
enum
{
MARKER_DIRAC = 0,
- MARKER_HEAVISIDE = 1
+ MARKER_HEAVISIDE = 1,
+ MARKER_NONE = 2
};
class FixPropertyAtomTracer : public FixPropertyAtom {
diff --git a/src/fix_property_global.cpp b/src/fix_property_global.cpp
index d106171..4edb8c0 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/fix_scalar_transport_equation.cpp b/src/fix_scalar_transport_equation.cpp
index 3cb2b91..37c7861 100644
--- a/src/fix_scalar_transport_equation.cpp
+++ b/src/fix_scalar_transport_equation.cpp
@@ -159,7 +159,7 @@ void FixScalarTransportEquation::post_create()
fixarg[5]="yes";
fixarg[6]="yes";
fixarg[7]="no";
- sprintf(fixarg[8],"%f",quantity_0);
+ sprintf(fixarg[8],"%e",quantity_0);
modify->add_fix(9,fixarg);
fix_quantity=static_cast<FixPropertyAtom*>(modify->find_fix_property(quantity_name,"property/atom","scalar",0,0,style));
}
diff --git a/src/general_container.h b/src/general_container.h
index ad17940..65d1e2f 100644
--- a/src/general_container.h
+++ b/src/general_container.h
@@ -132,7 +132,7 @@ namespace LAMMPS_NS
protected:
- GeneralContainer();
+ GeneralContainer(char *_id);
GeneralContainer(char *_id, char *_comm, char *_ref, char *_restart, int _scalePower = 1);
GeneralContainer(GeneralContainer<T,NUM_VEC,LEN_VEC> const &orig);
virtual ~GeneralContainer();
diff --git a/src/general_container_I.h b/src/general_container_I.h
index 4562bef..6eb36b8 100644
--- a/src/general_container_I.h
+++ b/src/general_container_I.h
@@ -33,8 +33,8 @@
------------------------------------------------------------------------- */
template<typename T, int NUM_VEC, int LEN_VEC>
- GeneralContainer<T,NUM_VEC,LEN_VEC>::GeneralContainer()
- : ContainerBase(),
+ GeneralContainer<T,NUM_VEC,LEN_VEC>::GeneralContainer(char *_id)
+ : ContainerBase(_id),
numElem_(0),
maxElem_(GROW)
{
@@ -604,12 +604,13 @@
template<typename T, int NUM_VEC, int LEN_VEC>
int GeneralContainer<T,NUM_VEC,LEN_VEC>::elemBufSize(int operation,bool scale,bool translate,bool rotate)
{
+
if(!this->decidePackUnpackOperation(operation,scale,translate,rotate))
return 0;
if(!this->decideCommOperation(operation))
return 0;
-
+
return (NUM_VEC*LEN_VEC);
}
diff --git a/src/lammps.cpp b/src/lammps.cpp
index 7d95660..e8500e2 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 3104dbe..9fb6adc 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 a5cda41..d9d619c 100644
--- a/src/modify.cpp
+++ b/src/modify.cpp
@@ -729,10 +729,14 @@ void Modify::add_fix(int narg, char **arg, char *suffix)
fmask[ifix] = fix[ifix]->setmask();
if (newflag) nfix++;
+ fix[ifix]->post_create_pre_restart();
+
// check if Fix is in restart_global list
// 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]);
@@ -744,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 c4852e7..f08099a 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 b03a086..4b014c6 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 67a4040..5d3e0e2 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 e12d152..d995e3c 100644
--- a/src/multi_node_mesh_I.h
+++ b/src/multi_node_mesh_I.h
@@ -35,6 +35,10 @@
template<int NUM_NODES>
MultiNodeMesh<NUM_NODES>::MultiNodeMesh(LAMMPS *lmp)
: AbstractMesh(lmp),
+ node_("node"),
+ nodesLastRe_("nodesLastRe"),
+ center_("center"),
+ rBound_("rBound"),
node_orig_(0),
precision_(EPSILON_PRECISION),
autoRemoveDuplicates_(false),
@@ -224,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;
}
@@ -245,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;
}
@@ -308,7 +326,7 @@
if(node_orig_)
error->one(FLERR,"Illegal situation in MultiNodeMesh<NUM_NODES>::registerMove");
- node_orig_ = new MultiVectorContainer<double,NUM_NODES,3>;
+ node_orig_ = new MultiVectorContainer<double,NUM_NODES,3>("node_orig");
for(int i = 0; i < nall; i++)
{
for(int j = 0; j < NUM_NODES; j++)
diff --git a/src/multi_node_mesh_parallel_buffer_I.h b/src/multi_node_mesh_parallel_buffer_I.h
index a0e04a4..3c84f8e 100644
--- a/src/multi_node_mesh_parallel_buffer_I.h
+++ b/src/multi_node_mesh_parallel_buffer_I.h
@@ -253,7 +253,7 @@
int nsend = 0;
- if(operation == OPERATION_COMM_EXCHANGE || operation == OPERATION_COMM_BORDERS)
+ if(OPERATION_COMM_EXCHANGE == operation || OPERATION_COMM_BORDERS == operation)
{
nsend += MultiNodeMesh<NUM_NODES>::center_.pushElemListToBuffer(n,list,&(buf[nsend]),operation);
@@ -264,7 +264,7 @@
return nsend;
}
- if(operation == OPERATION_COMM_FORWARD)
+ if(OPERATION_COMM_FORWARD == operation)
{
return nsend;
@@ -285,7 +285,7 @@
{
int nrecv = 0;
- if(operation == OPERATION_COMM_EXCHANGE || operation == OPERATION_COMM_BORDERS)
+ if(OPERATION_COMM_EXCHANGE == operation || OPERATION_COMM_BORDERS == operation)
{
nrecv += MultiNodeMesh<NUM_NODES>::center_.popElemListFromBuffer(first,n,&(buf[nrecv]),operation);
nrecv += MultiNodeMesh<NUM_NODES>::node_.popElemListFromBuffer(first,n,&(buf[nrecv]),operation);
@@ -295,7 +295,7 @@
return nrecv;
}
- if(operation == OPERATION_COMM_FORWARD)
+ if(OPERATION_COMM_FORWARD == operation)
{
// nrecv += MultiNodeMesh<NUM_NODES>::node_.popListFromBuffer(first,n,&(buf[nrecv]),operation);
@@ -317,7 +317,7 @@
{
int nsend = 0;
- if(operation == OPERATION_COMM_REVERSE)
+ if(OPERATION_COMM_REVERSE == operation)
{
return nsend;
@@ -338,7 +338,7 @@
{
int nrecv = 0;
- if(operation == OPERATION_COMM_REVERSE)
+ if(OPERATION_COMM_REVERSE == operation)
{
return nrecv;
@@ -360,13 +360,13 @@
{
int size_buf = 0;
- if(operation == OPERATION_RESTART)
+ if(OPERATION_RESTART == operation)
{
size_buf += MultiNodeMesh<NUM_NODES>::node_.elemBufSize();
return size_buf;
}
- if(operation == OPERATION_COMM_EXCHANGE || operation == OPERATION_COMM_BORDERS)
+ if(OPERATION_COMM_EXCHANGE == operation || OPERATION_COMM_BORDERS == operation)
{
size_buf += MultiNodeMesh<NUM_NODES>::center_.elemBufSize();
size_buf += MultiNodeMesh<NUM_NODES>::node_.elemBufSize();
@@ -376,13 +376,13 @@
return size_buf;
}
- if(operation == OPERATION_COMM_FORWARD)
+ if(OPERATION_COMM_FORWARD == operation)
{
return size_buf;
}
- if(operation == OPERATION_COMM_REVERSE)
+ if(OPERATION_COMM_REVERSE == operation)
{
return size_buf;
@@ -403,14 +403,14 @@
{
int nsend = 0;
- if(operation == OPERATION_RESTART)
+ if(OPERATION_RESTART == operation)
{
nsend += MultiNodeMesh<NUM_NODES>::node_.pushElemToBuffer(i,&(buf[nsend]),operation);
return nsend;
}
- if(operation == OPERATION_COMM_EXCHANGE || operation == OPERATION_COMM_BORDERS)
+ if(OPERATION_COMM_EXCHANGE == operation || OPERATION_COMM_BORDERS == operation)
{
nsend += MultiNodeMesh<NUM_NODES>::center_.pushElemToBuffer(i,&(buf[nsend]),operation);
@@ -436,9 +436,9 @@
{
int nrecv = 0;
- if(operation == OPERATION_RESTART)
+ if(OPERATION_RESTART == operation)
{
- MultiVectorContainer<double,NUM_NODES,3> nodeTmp;
+ MultiVectorContainer<double,NUM_NODES,3> nodeTmp("nodeTmp");
nrecv += nodeTmp.popElemFromBuffer(&(buf[nrecv]),operation);
this->addElement(nodeTmp.begin()[0],-1);
@@ -448,7 +448,7 @@
return nrecv;
}
- if(operation == OPERATION_COMM_EXCHANGE || operation == OPERATION_COMM_BORDERS)
+ if(OPERATION_COMM_EXCHANGE == operation || OPERATION_COMM_BORDERS == operation)
{
nrecv += MultiNodeMesh<NUM_NODES>::center_.popElemFromBuffer(&(buf[nrecv]),operation);
nrecv += MultiNodeMesh<NUM_NODES>::node_.popElemFromBuffer(&(buf[nrecv]),operation);
diff --git a/src/multi_vector_container.h b/src/multi_vector_container.h
index 866d475..11b7338 100644
--- a/src/multi_vector_container.h
+++ b/src/multi_vector_container.h
@@ -35,7 +35,7 @@ namespace LAMMPS_NS{
class MultiVectorContainer : public GeneralContainer <T, NUM_VEC, LEN_VEC>
{
public:
- MultiVectorContainer();
+ MultiVectorContainer(char *_id);
MultiVectorContainer(char *_id, char *_comm, char *_ref, char *_restart, int _scalePower = 1);
MultiVectorContainer(MultiVectorContainer<T,NUM_VEC,LEN_VEC> const &orig);
virtual ~MultiVectorContainer();
@@ -46,8 +46,8 @@ namespace LAMMPS_NS{
------------------------------------------------------------------------- */
template<typename T, int NUM_VEC, int LEN_VEC>
- MultiVectorContainer<T,NUM_VEC,LEN_VEC>::MultiVectorContainer()
- : GeneralContainer<T,NUM_VEC,LEN_VEC>()
+ MultiVectorContainer<T,NUM_VEC,LEN_VEC>::MultiVectorContainer(char *_id)
+ : GeneralContainer<T,NUM_VEC,LEN_VEC>(_id)
{
}
diff --git a/src/neigh_gran.cpp b/src/neigh_gran.cpp
index 5c08b21..8de54ae 100644
--- a/src/neigh_gran.cpp
+++ b/src/neigh_gran.cpp
@@ -73,6 +73,8 @@ void Neighbor::granular_nsq_no_newton(NeighList *list)
int **firstneigh = list->firstneigh;
int **pages = list->pages;
+ double contactHistoryDistanceFactorSqr = contactHistoryDistanceFactor*contactHistoryDistanceFactor;
+
FixContactHistory *fix_history = list->fix_history;
if (fix_history) {
npartner = fix_history->npartner;
@@ -134,8 +136,8 @@ void Neighbor::granular_nsq_no_newton(NeighList *list)
neighptr[n] = j;
if (fix_history) {
- if (rsq < radsum*radsum*contactHistoryDistanceFactor)
- {
+ if (rsq < radsum*radsum*contactHistoryDistanceFactorSqr)
+ {
for (m = 0; m < npartner[i]; m++)
if (partner[i][m] == tag[j]) break;
if (m < npartner[i]) {
@@ -300,6 +302,8 @@ void Neighbor::granular_bin_no_newton(NeighList *list)
double **pages_shear;
int dnum;
+ double contactHistoryDistanceFactorSqr = contactHistoryDistanceFactor*contactHistoryDistanceFactor;
+
// bin local & ghost atoms
bin_atoms();
@@ -383,12 +387,12 @@ 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;
if (fix_history) {
- if (rsq < radsum*radsum*contactHistoryDistanceFactor)
+ if (rsq < radsum*radsum*contactHistoryDistanceFactorSqr)
{
for (m = 0; m < npartner[i]; m++)
if (partner[i][m] == tag[j]) break;
diff --git a/src/neighbor.cpp b/src/neighbor.cpp
index 6ed0325..8fd321c 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;
}
@@ -1682,6 +1685,11 @@ void Neighbor::modify_params(int narg, char **arg)
delay = atoi(arg[iarg+1]);
if (delay < 0) error->all(FLERR,"Illegal neigh_modify command");
iarg += 2;
+ } else if (strcmp(arg[iarg],"contactHistoryDistanceFactor") == 0) {
+ if (iarg+2 > narg) error->all(FLERR,"Illegal neigh_modify command");
+ contactHistoryDistanceFactor = atof(arg[iarg+1]);
+ if (contactHistoryDistanceFactor < 1.0) error->all(FLERR,"Illegal neigh_modify command. Please set contactHistoryDistanceFactor value >=1");
+ iarg +=2;
} else if (strcmp(arg[iarg],"check") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal neigh_modify command");
if (strcmp(arg[iarg+1],"yes") == 0) dist_check = 1;
diff --git a/src/pair_gran.cpp b/src/pair_gran.cpp
index 79f3cf4..6e2c7eb 100644
--- a/src/pair_gran.cpp
+++ b/src/pair_gran.cpp
@@ -192,6 +192,9 @@ void PairGran::init_style()
// error and warning checks
+ 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)
error->all(FLERR,"Cannot use a non-consistent unit system with pair gran. Please use si,cgs or lj.");
diff --git a/src/pair_gran_hooke.cpp b/src/pair_gran_hooke.cpp
index 4594f00..249a536 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 8cfc5b1..f5463dd 100644
--- a/src/pair_gran_hooke_history.cpp
+++ b/src/pair_gran_hooke_history.cpp
@@ -663,7 +663,7 @@ void PairGranHookeHistory::init_granular()
if(rollingflag)
coeffRollFrict1=static_cast<FixPropertyGlobal*>(modify->find_fix_property("coefficientRollingFriction","property/global","peratomtypepair",max_type,max_type,force->pair_style));
- if(rollingflag == 2 || rollingflag == 3) // epsd model
+ if(rollingflag == 2) // damping for original epsd model only
coeffRollVisc1=static_cast<FixPropertyGlobal*>(modify->find_fix_property("coefficientRollingViscousDamping","property/global","peratomtypepair",max_type,max_type,force->pair_style));
if(viscousflag)
{
@@ -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++)
@@ -731,7 +732,7 @@ void PairGranHookeHistory::init_granular()
coeffFrict[i][j] = coeffFrict1->compute_array(i-1,j-1);
if(rollingflag) coeffRollFrict[i][j] = coeffRollFrict1->compute_array(i-1,j-1);
- if(rollingflag == 2 || rollingflag == 3) coeffRollVisc[i][j] = coeffRollVisc1->compute_array(i-1,j-1);
+ if(rollingflag == 2) coeffRollVisc[i][j] = coeffRollVisc1->compute_array(i-1,j-1);
if(cohesionflag) cohEnergyDens[i][j] = cohEnergyDens1->compute_array(i-1,j-1);
//omitting veff here
@@ -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/pair_hybrid.cpp b/src/pair_hybrid.cpp
index a31f864..5ae4f2f 100644
--- a/src/pair_hybrid.cpp
+++ b/src/pair_hybrid.cpp
@@ -253,6 +253,16 @@ void PairHybrid::settings(int narg, char **arg)
// set comm_forward, comm_reverse, comm_reverse_off to max of any sub-style
+ flags();
+}
+
+/* ----------------------------------------------------------------------
+ set top-level pair flags from sub-style flags - new from patch 19Jan2013
+------------------------------------------------------------------------- */
+
+void PairHybrid::flags() //- new from patch 19Jan2013
+{
+ int m;
for (m = 0; m < nstyles; m++) {
if (styles[m]) comm_forward = MAX(comm_forward,styles[m]->comm_forward);
if (styles[m]) comm_reverse = MAX(comm_reverse,styles[m]->comm_reverse);
@@ -602,6 +612,7 @@ void PairHybrid::read_restart(FILE *fp)
styles = new Pair*[nstyles];
keywords = new char*[nstyles];
+ multiple = new int[nstyles]; // - new from patch 17May2012
// each sub-style is created via new_pair()
// each reads its settings, but no coeff info
@@ -615,6 +626,18 @@ void PairHybrid::read_restart(FILE *fp)
styles[m] = force->new_pair(keywords[m],lmp->suffix,dummy);
styles[m]->read_restart_settings(fp);
}
+ for (int i = 0; i < nstyles; i++) {
+ int count = 0;
+ for (int j = 0; j < nstyles; j++) {
+ if (strcmp(keywords[j],keywords[i]) == 0) count++;
+ if (j == i) multiple[i] = count;
+ }
+ if (count == 1) multiple[i] = 0;
+ } //- new from patch 17May2012
+
+ // set pair flags from sub-style flags // - new from patch 19Jan2013
+
+ flags();
}
/* ----------------------------------------------------------------------
diff --git a/src/pair_hybrid.h b/src/pair_hybrid.h
index d196498..913bb7c 100644
--- a/src/pair_hybrid.h
+++ b/src/pair_hybrid.h
@@ -61,6 +61,7 @@ class PairHybrid : public Pair {
int nallstyles;
void allocate();
+ void flags();
virtual void modify_requests();
void build_styles();
int known_style(char *);
diff --git a/src/pair_sph.cpp b/src/pair_sph.cpp
index cba9a7c..e788ad3 100644
--- a/src/pair_sph.cpp
+++ b/src/pair_sph.cpp
@@ -68,7 +68,27 @@ PairSph::PairSph(LAMMPS *lmp) : Pair(lmp)
sl = NULL;
slComType = NULL;
+ fix_fgradP_ = NULL;
+
mass_type = atom->avec->mass_type; // get flag for mass per type
+
+ char **fixarg;
+ fixarg=new char*[11];
+ for (int kk=0;kk<11;kk++) fixarg[kk]=new char[30];
+
+ strcpy(fixarg[0],"fgradP");
+ fixarg[1]="all";
+ fixarg[2]="property/atom";
+ strcpy(fixarg[3],"fgradP");
+ fixarg[4]="vector";
+ fixarg[5]="yes";
+ fixarg[6]="yes";
+ fixarg[7]="yes";
+ fixarg[8]="0.";
+ fixarg[9]="0.";
+ fixarg[10]="0.";
+ fix_fgradP_ = modify->add_fix_property_atom(11,fixarg,"PairSph");
+ delete []fixarg;
}
/* ---------------------------------------------------------------------- */
@@ -87,6 +107,8 @@ PairSph::~PairSph()
if(kernel_style) delete []kernel_style;
if(fppaSl) modify->delete_fix("sl");
// if(fppaSlType) modify->delete_fix("sl");
+
+ if(fix_fgradP_) modify->delete_fix(fix_fgradP_->id);
}
/* ----------------------------------------------------------------------
diff --git a/src/pair_sph.h b/src/pair_sph.h
index 3a7a60d..b4d8a8e 100644
--- a/src/pair_sph.h
+++ b/src/pair_sph.h
@@ -91,6 +91,9 @@ class PairSph : public Pair {
int pairStyle_;
double viscosity_;
+ // storage for force part caused by pressure gradient (grad P / rho):
+ class FixPropertyAtom* fix_fgradP_;
+ double **fgradP_;
};
}
diff --git a/src/pair_sph_artvisc_tenscorr.cpp b/src/pair_sph_artvisc_tenscorr.cpp
index f24b3ba..8c71999 100644
--- a/src/pair_sph_artvisc_tenscorr.cpp
+++ b/src/pair_sph_artvisc_tenscorr.cpp
@@ -472,6 +472,10 @@ void PairSphArtviscTenscorr::compute_eval(int eflag, int vflag)
// get distance and normalized distance
r = sqrt(rsq);
+ if (r == 0.) {
+ printf("Particle %i and %i are at same position (%f, %f, %f)",i,j,xtmp,ytmp,ztmp);
+ error->one(FLERR,"Zero distance between SPH particles!");
+ }
rinv = 1./r;
s = r * slComInv;
diff --git a/src/pointers.h b/src/pointers.h
index 981b5cc..6fce2ed 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/primitive_wall.h b/src/primitive_wall.h
index 6955446..cb7d55b 100644
--- a/src/primitive_wall.h
+++ b/src/primitive_wall.h
@@ -72,7 +72,7 @@ namespace LAMMPS_NS
*/
PrimitiveWall::PrimitiveWall(LAMMPS *lmp,PRIMITIVE_WALL_DEFINITIONS::WallType wType_, int nParam_, double *param_)
- : Pointers(lmp), wType(wType_), nParam(nParam_)
+ : Pointers(lmp), neighlist("neighlist"), wType(wType_), nParam(nParam_)
{
param = new double[nParam];
for(int i=0;i<nParam;i++)
diff --git a/src/read_data.cpp b/src/read_data.cpp
index 4075e90..4b6eaec 100644
--- a/src/read_data.cpp
+++ b/src/read_data.cpp
@@ -97,7 +97,7 @@ ReadData::~ReadData()
void ReadData::command(int narg, char **arg)
{
if (narg != 1 && narg != 2) error->all(FLERR,"Illegal read_data command");
-
+
if (narg == 2 && strcmp(arg[1],"add") == 0) add_to_existing = 1;
if (domain->box_exist && !add_to_existing)
@@ -150,10 +150,9 @@ void ReadData::command(int narg, char **arg)
}
}
- header(1,atom->molecular?0:1);
+ header(1,(atom->molecular?0:1) || (0 < comm->me));
if (!domain->box_exist) domain->box_exist = 1;
-
// problem setup using info from header
if(!add_to_existing)
@@ -527,8 +526,7 @@ void ReadData::header(int flag, int add)
atom->nangles < 0 || atom->nangles > MAXBIGINT ||
atom->ndihedrals < 0 || atom->ndihedrals > MAXBIGINT ||
atom->nimpropers < 0 || atom->nimpropers > MAXBIGINT) {
- if (flag == 0) error->one(FLERR,"System in data file is too big");
- else error->all(FLERR,"System in data file is too big");
+ if (me == 0) error->one(FLERR,"System in data file is too big"); //- new from patch 17May2012
}
// check that exiting string is a valid section keyword
@@ -536,34 +534,34 @@ void ReadData::header(int flag, int add)
parse_keyword(1,flag);
for (n = 0; n < NSECTIONS; n++)
if (strcmp(keyword,section_keywords[n]) == 0) break;
- if (n == NSECTIONS) {
+ if (n == NSECTIONS && me == 0) { //- new from patch 17May2012
char str[128];
sprintf(str,"Unknown identifier in data file: %s",keyword);
- error->all(FLERR,str);
+ error->one(FLERR,str); //- new from patch 17May2012
}
// error check on consistency of header values
if ((atom->nbonds || atom->nbondtypes) &&
- atom->avec->bonds_allow == 0)
+ atom->avec->bonds_allow == 0 && me == 0) //- new from patch 17May2012
error->one(FLERR,"No bonds allowed with this atom style");
if ((atom->nangles || atom->nangletypes) &&
- atom->avec->angles_allow == 0)
+ atom->avec->angles_allow == 0 && me == 0) //- new from patch 17May2012
error->one(FLERR,"No angles allowed with this atom style");
if ((atom->ndihedrals || atom->ndihedraltypes) &&
- atom->avec->dihedrals_allow == 0)
+ atom->avec->dihedrals_allow == 0 && me == 0) //- new from patch 17May2012
error->one(FLERR,"No dihedrals allowed with this atom style");
if ((atom->nimpropers || atom->nimpropertypes) &&
- atom->avec->impropers_allow == 0)
+ atom->avec->impropers_allow == 0 && me == 0) //- new from patch 17May2012
error->one(FLERR,"No impropers allowed with this atom style");
- if (atom->nbonds > 0 && atom->nbondtypes <= 0)
+ if (atom->nbonds > 0 && atom->nbondtypes <= 0 && me == 0) //- new from patch 17May2012
error->one(FLERR,"Bonds defined but no bond types");
- if (atom->nangles > 0 && atom->nangletypes <= 0)
+ if (atom->nangles > 0 && atom->nangletypes <= 0 && me == 0) //- new from patch 17May2012
error->one(FLERR,"Angles defined but no angle types");
- if (atom->ndihedrals > 0 && atom->ndihedraltypes <= 0)
+ if (atom->ndihedrals > 0 && atom->ndihedraltypes <= 0 && me == 0) //- new from patch 17May2012
error->one(FLERR,"Dihedrals defined but no dihedral types");
- if (atom->nimpropers > 0 && atom->nimpropertypes <= 0)
+ if (atom->nimpropers > 0 && atom->nimpropertypes <= 0 && me == 0) //- new from patch 17May2012
error->one(FLERR,"Impropers defined but no improper types");
}
@@ -1195,107 +1193,107 @@ void ReadData::scan(int &bond_per_atom, int &angle_per_atom,
memory->create(count,cmax,"read_data:count");
while (strlen(keyword)) {
-
+
if (strcmp(keyword,"Masses") == 0) skip_lines(atom->ntypes);
- else if (strcmp(keyword,"Atoms") == 0) skip_lines(natoms);
- else if (strcmp(keyword,"Velocities") == 0) skip_lines(natoms);
+ else if (strcmp(keyword,"Atoms") == 0) skip_lines(add_to_existing?natoms_add:natoms);
+ else if (strcmp(keyword,"Velocities") == 0) skip_lines(add_to_existing?natoms_add:natoms);
else if (strcmp(keyword,"Ellipsoids") == 0) {
if (!avec_ellipsoid)
- error->all(FLERR,"Invalid data file section: Ellipsoids");
+ error->one(FLERR,"Invalid data file section: Ellipsoids");
ellipsoid_flag = 1;
skip_lines(nellipsoids);
} else if (strcmp(keyword,"Lines") == 0) {
- if (!avec_line) error->all(FLERR,"Invalid data file section: Lines");
+ if (!avec_line) error->one(FLERR,"Invalid data file section: Lines");
line_flag = 1;
skip_lines(nlines);
} else if (strcmp(keyword,"Triangles") == 0) {
- if (!avec_tri) error->all(FLERR,"Invalid data file section: Triangles");
+ if (!avec_tri) error->one(FLERR,"Invalid data file section: Triangles");
tri_flag = 1;
skip_lines(ntris);
} else if (strcmp(keyword,"Pair Coeffs") == 0) {
if (force->pair == NULL)
- error->all(FLERR,"Must define pair_style before Pair Coeffs");
+ error->one(FLERR,"Must define pair_style before Pair Coeffs");
skip_lines(atom->ntypes);
} else if (strcmp(keyword,"Bond Coeffs") == 0) {
if (atom->avec->bonds_allow == 0)
- error->all(FLERR,"Invalid data file section: Bond Coeffs");
+ error->one(FLERR,"Invalid data file section: Bond Coeffs");
if (force->bond == NULL)
- error->all(FLERR,"Must define bond_style before Bond Coeffs");
+ error->one(FLERR,"Must define bond_style before Bond Coeffs");
skip_lines(atom->nbondtypes);
} else if (strcmp(keyword,"Angle Coeffs") == 0) {
if (atom->avec->angles_allow == 0)
- error->all(FLERR,"Invalid data file section: Angle Coeffs");
+ error->one(FLERR,"Invalid data file section: Angle Coeffs");
if (force->angle == NULL)
- error->all(FLERR,"Must define angle_style before Angle Coeffs");
+ error->one(FLERR,"Must define angle_style before Angle Coeffs");
skip_lines(atom->nangletypes);
} else if (strcmp(keyword,"Dihedral Coeffs") == 0) {
skip_lines(atom->ndihedraltypes);
if (atom->avec->dihedrals_allow == 0)
- error->all(FLERR,"Invalid data file section: Dihedral Coeffs");
+ error->one(FLERR,"Invalid data file section: Dihedral Coeffs");
if (force->dihedral == NULL)
- error->all(FLERR,"Must define dihedral_style before Dihedral Coeffs");
+ error->one(FLERR,"Must define dihedral_style before Dihedral Coeffs");
} else if (strcmp(keyword,"Improper Coeffs") == 0) {
if (atom->avec->impropers_allow == 0)
- error->all(FLERR,"Invalid data file section: Improper Coeffs");
+ error->one(FLERR,"Invalid data file section: Improper Coeffs");
if (force->improper == NULL)
- error->all(FLERR,"Must define improper_style before Improper Coeffs");
+ error->one(FLERR,"Must define improper_style before Improper Coeffs");
skip_lines(atom->nimpropertypes);
} else if (strcmp(keyword,"BondBond Coeffs") == 0) {
if (atom->avec->angles_allow == 0)
- error->all(FLERR,"Invalid data file section: BondBond Coeffs");
+ error->one(FLERR,"Invalid data file section: BondBond Coeffs");
if (force->angle == NULL)
- error->all(FLERR,"Must define angle_style before BondBond Coeffs");
+ error->one(FLERR,"Must define angle_style before BondBond Coeffs");
skip_lines(atom->nangletypes);
} else if (strcmp(keyword,"BondAngle Coeffs") == 0) {
if (atom->avec->angles_allow == 0)
- error->all(FLERR,"Invalid data file section: BondAngle Coeffs");
+ error->one(FLERR,"Invalid data file section: BondAngle Coeffs");
if (force->angle == NULL)
- error->all(FLERR,"Must define angle_style before BondAngle Coeffs");
+ error->one(FLERR,"Must define angle_style before BondAngle Coeffs");
skip_lines(atom->nangletypes);
} else if (strcmp(keyword,"MiddleBondTorsion Coeffs") == 0) {
if (atom->avec->dihedrals_allow == 0)
- error->all(FLERR,"Invalid data file section: MiddleBondTorsion Coeffs");
+ error->one(FLERR,"Invalid data file section: MiddleBondTorsion Coeffs");
if (force->dihedral == NULL)
- error->all(FLERR,
+ error->one(FLERR,
"Must define dihedral_style before "
"MiddleBondTorsion Coeffs");
skip_lines(atom->ndihedraltypes);
} else if (strcmp(keyword,"EndBondTorsion Coeffs") == 0) {
if (atom->avec->dihedrals_allow == 0)
- error->all(FLERR,"Invalid data file section: EndBondTorsion Coeffs");
+ error->one(FLERR,"Invalid data file section: EndBondTorsion Coeffs");
if (force->dihedral == NULL)
- error->all(FLERR,
+ error->one(FLERR,
"Must define dihedral_style before EndBondTorsion Coeffs");
skip_lines(atom->ndihedraltypes);
} else if (strcmp(keyword,"AngleTorsion Coeffs") == 0) {
if (atom->avec->dihedrals_allow == 0)
- error->all(FLERR,"Invalid data file section: AngleTorsion Coeffs");
+ error->one(FLERR,"Invalid data file section: AngleTorsion Coeffs");
if (force->dihedral == NULL)
- error->all(FLERR,
+ error->one(FLERR,
"Must define dihedral_style before AngleTorsion Coeffs");
skip_lines(atom->ndihedraltypes);
} else if (strcmp(keyword,"AngleAngleTorsion Coeffs") == 0) {
if (atom->avec->dihedrals_allow == 0)
- error->all(FLERR,"Invalid data file section: AngleAngleTorsion Coeffs");
+ error->one(FLERR,"Invalid data file section: AngleAngleTorsion Coeffs");
if (force->dihedral == NULL)
- error->all(FLERR,
+ error->one(FLERR,
"Must define dihedral_style before "
"AngleAngleTorsion Coeffs");
skip_lines(atom->ndihedraltypes);
} else if (strcmp(keyword,"BondBond13 Coeffs") == 0) {
if (atom->avec->dihedrals_allow == 0)
- error->all(FLERR,"Invalid data file section: BondBond13 Coeffs");
+ error->one(FLERR,"Invalid data file section: BondBond13 Coeffs");
if (force->dihedral == NULL)
- error->all(FLERR,"Must define dihedral_style before BondBond13 Coeffs");
+ error->one(FLERR,"Must define dihedral_style before BondBond13 Coeffs");
skip_lines(atom->ndihedraltypes);
} else if (strcmp(keyword,"AngleAngle Coeffs") == 0) {
if (atom->avec->impropers_allow == 0)
- error->all(FLERR,"Invalid data file section: AngleAngle Coeffs");
+ error->one(FLERR,"Invalid data file section: AngleAngle Coeffs");
if (force->improper == NULL)
- error->all(FLERR,"Must define improper_style before AngleAngle Coeffs");
+ error->one(FLERR,"Must define improper_style before AngleAngle Coeffs");
skip_lines(atom->nimpropertypes);
} else if (strcmp(keyword,"Bonds") == 0) {
diff --git a/src/region.cpp b/src/region.cpp
index 45ec3f1..0e7beb8 100644
--- a/src/region.cpp
+++ b/src/region.cpp
@@ -586,10 +586,10 @@ void Region::volume_mc(int n_test,bool cutflag,double cut,double &vol_global,dou
}
else
{
- if(match(pos[0],pos[1],pos[2]) && !match_cut(pos,cut) )
+ if(match(pos[0],pos[1],pos[2]))
{
n_in_global++;
- if(domain->is_in_subdomain(pos))
+ if(domain->is_in_subdomain(pos) && !match_cut(pos,cut) )
n_in_local++;
}
}
@@ -608,10 +608,9 @@ void Region::volume_mc(int n_test,bool cutflag,double cut,double &vol_global,dou
vol_global = static_cast<double>(n_in_global_all)/static_cast<double>(n_test*comm->nprocs) * vol_bbox;
vol_local = static_cast<double>(n_in_local )/static_cast<double>(n_test) * vol_bbox;
- // sum of local volumes will not be equal to global volume because of
- // different random generator states - correct this now
MPI_Sum_Scalar(vol_local,vol_local_all,world);
vol_local *= (vol_global/vol_local_all);
+
}
/* ---------------------------------------------------------------------- */
diff --git a/src/region_wedge.cpp b/src/region_wedge.cpp
index 88ab742..0c7e0af 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 3164422..0d0749b 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/scalar_container.h b/src/scalar_container.h
index b141d24..78f3d69 100644
--- a/src/scalar_container.h
+++ b/src/scalar_container.h
@@ -38,7 +38,7 @@ namespace LAMMPS_NS
class ScalarContainer : public GeneralContainer <T, 1, 1>
{
public:
- ScalarContainer();
+ ScalarContainer(char *_id);
ScalarContainer(char *_id, char *_comm, char *_ref, char *_restart, int _scalePower = 1);
ScalarContainer(ScalarContainer<T> const &orig);
virtual ~ScalarContainer();
@@ -60,8 +60,8 @@ namespace LAMMPS_NS
------------------------------------------------------------------------- */
template<typename T>
- ScalarContainer<T>::ScalarContainer()
- : GeneralContainer<T,1,1>()
+ ScalarContainer<T>::ScalarContainer(char *_id)
+ : GeneralContainer<T,1,1>(_id)
{
}
diff --git a/src/set.cpp b/src/set.cpp
index 9bcef1f..2001070 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 de8579e..2358369 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 a175268..7085f6c 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_container.h b/src/vector_container.h
index 3948137..7927e91 100644
--- a/src/vector_container.h
+++ b/src/vector_container.h
@@ -37,7 +37,7 @@ namespace LAMMPS_NS
class VectorContainer : public GeneralContainer <T, 1, LEN_VEC>
{
public:
- VectorContainer();
+ VectorContainer(char *_id);
VectorContainer(char *_id, char *_comm, char *_ref, char *_restart, int _scalePower = 1);
VectorContainer(VectorContainer<T,LEN_VEC> const &orig);
virtual ~VectorContainer();
@@ -62,8 +62,8 @@ namespace LAMMPS_NS
------------------------------------------------------------------------- */
template<typename T, int LEN_VEC>
- VectorContainer<T,LEN_VEC>::VectorContainer()
- : GeneralContainer<T,1,LEN_VEC>()
+ VectorContainer<T,LEN_VEC>::VectorContainer(char *_id)
+ : GeneralContainer<T,1,LEN_VEC>(_id)
{
}
diff --git a/src/vector_liggghts.h b/src/vector_liggghts.h
index a8f17b5..2ca74f6 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 758600f..dbaa011 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 8b9ba5f..53eb557 100644
--- a/src/version_liggghts.h
+++ b/src/version_liggghts.h
@@ -1 +1 @@
-#define LIGGGHTS_VERSION "LIGGGHTS-PUBLIC 2.3.6, compiled 2013-07-02-17:35:42 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 e75da3e..bc4abe8 100644
--- a/src/version_liggghts.txt
+++ b/src/version_liggghts.txt
@@ -1 +1 @@
-2.3.6
+2.3.8
diff --git a/src/volume_mesh_I.h b/src/volume_mesh_I.h
index fe4bafd..20fae8b 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))
{
--
Alioth's /usr/local/bin/git-commit-notice on /srv/git.debian.org/git/debian-science/packages/liggghts.git
More information about the debian-science-commits
mailing list