1 /* ----------------------------------------------------------------------
2     This is the
3 
4     ██╗     ██╗ ██████╗  ██████╗  ██████╗ ██╗  ██╗████████╗███████╗
5     ██║     ██║██╔════╝ ██╔════╝ ██╔════╝ ██║  ██║╚══██╔══╝██╔════╝
6     ██║     ██║██║  ███╗██║  ███╗██║  ███╗███████║   ██║   ███████╗
7     ██║     ██║██║   ██║██║   ██║██║   ██║██╔══██║   ██║   ╚════██║
8     ███████╗██║╚██████╔╝╚██████╔╝╚██████╔╝██║  ██║   ██║   ███████║
9     ╚══════╝╚═╝ ╚═════╝  ╚═════╝  ╚═════╝ ╚═╝  ╚═╝   ╚═╝   ╚══════╝®
10 
11     DEM simulation engine, released by
12     DCS Computing Gmbh, Linz, Austria
13     http://www.dcs-computing.com, office@dcs-computing.com
14 
15     LIGGGHTS® is part of CFDEM®project:
16     http://www.liggghts.com | http://www.cfdem.com
17 
18     Core developer and main author:
19     Christoph Kloss, christoph.kloss@dcs-computing.com
20 
21     LIGGGHTS® is open-source, distributed under the terms of the GNU Public
22     License, version 2 or later. It is distributed in the hope that it will
23     be useful, but WITHOUT ANY WARRANTY; without even the implied warranty
24     of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. You should have
25     received a copy of the GNU General Public License along with LIGGGHTS®.
26     If not, see http://www.gnu.org/licenses . See also top-level README
27     and LICENSE files.
28 
29     LIGGGHTS® and CFDEM® are registered trade marks of DCS Computing GmbH,
30     the producer of the LIGGGHTS® software and the CFDEM®coupling software
31     See http://www.cfdem.com/terms-trademark-policy for details.
32 
33 -------------------------------------------------------------------------
34     Contributing author and copyright for this file:
35     (if not contributing author is listed, this file has been contributed
36     by the core developer)
37 
38     Copyright 2012-     DCS Computing GmbH, Linz
39     Copyright 2009-2012 JKU Linz
40 ------------------------------------------------------------------------- */
41 
42 #ifndef LMP_DOMAIN_I_H
43 #define LMP_DOMAIN_I_H
44 
45 /* ----------------------------------------------------------------------
46    check if coordinate in domain, subdomain or extended subdomain
47    need to test with <= and >= for domain, and < and >= for subdomain
48    inlined for performance
49 ------------------------------------------------------------------------- */
50 
is_in_domain(double * pos)51 inline int Domain::is_in_domain(double* pos)
52 {
53     if(is_wedge)
54         return is_in_domain_wedge(pos);
55 
56     if
57     (
58         pos[0] >= boxlo[0] && pos[0] <= boxhi[0] &&
59         pos[1] >= boxlo[1] && pos[1] <= boxhi[1] &&
60         pos[2] >= boxlo[2] && pos[2] <= boxhi[2]
61     )   return 1;
62     return 0;
63 }
64 
is_in_subdomain(double * pos)65 inline int Domain::is_in_subdomain(double* pos)
66 {
67     if(is_wedge)
68         return is_in_subdomain_wedge(pos);
69 
70     double checkhi[3];
71     double checklo[3];
72 
73     // If subdomain touches the lower or upper boundaries of the bounding box then
74     // add a small padding to avoid rounding errors
75     checkhi[0] = subhi[0] + (MathExtraLiggghts::compDouble(subhi[0], boxhi[0]) ? SMALL_DMBRDR : 0.0);
76     checkhi[1] = subhi[1] + (MathExtraLiggghts::compDouble(subhi[1], boxhi[1]) ? SMALL_DMBRDR : 0.0);
77     checkhi[2] = subhi[2] + (MathExtraLiggghts::compDouble(subhi[2], boxhi[2]) ? SMALL_DMBRDR : 0.0);
78     checklo[0] = sublo[0] - (MathExtraLiggghts::compDouble(sublo[0], boxlo[0]) ? SMALL_DMBRDR : 0.0);
79     checklo[1] = sublo[1] - (MathExtraLiggghts::compDouble(sublo[1], boxlo[1]) ? SMALL_DMBRDR : 0.0);
80     checklo[2] = sublo[2] - (MathExtraLiggghts::compDouble(sublo[2], boxlo[2]) ? SMALL_DMBRDR : 0.0);
81 
82     if ( pos[0] >= checklo[0] && pos[0] < checkhi[0] &&
83          pos[1] >= checklo[1] && pos[1] < checkhi[1] &&
84          pos[2] >= checklo[2] && pos[2] < checkhi[2])
85         return 1;
86     return 0;
87 }
88 
is_in_extended_subdomain(double * pos)89 inline int Domain::is_in_extended_subdomain(double* pos)
90 {
91     if(is_wedge)
92         return is_in_extended_subdomain_wedge(pos);
93 
94     // called on insertion
95     // yields true if particle would be in subdomain after box extension
96 
97     if (is_in_subdomain(pos))
98         return 1;
99     else if (dimension == 2)
100         error->all(FLERR,"Domain::is_in_extended_subdomain() not implemented for 2d");
101     else
102     {
103         bool flag = true;
104         for(int idim = 0; idim < 3; idim++)
105         {
106 
107             if (comm->procgrid[idim] == 1) {}
108             else if(comm->myloc[idim] == comm->procgrid[idim]-1)
109                 flag = flag && (pos[idim] >= sublo[idim]);
110             else if(comm->myloc[idim] == 0)
111                 flag = flag && (pos[idim] <= subhi[idim]);
112 
113             else
114                 flag = flag && (pos[idim] >= sublo[idim] && pos[idim] < subhi[idim]);
115         }
116         if(flag) return 1;
117         return 0;
118     }
119     return 0;
120 }
121 
122 /* ----------------------------------------------------------------------
123    check distance from borders of subbox
124 ------------------------------------------------------------------------- */
125 
dist_subbox_borders(double * pos)126 inline double Domain::dist_subbox_borders(double* pos)
127 {
128     if(is_wedge)
129         return dist_subbox_borders_wedge(pos);
130 
131     double deltalo[3], deltahi[3], dlo, dhi;
132 
133     vectorSubtract3D(sublo,pos,deltalo);
134     vectorSubtract3D(subhi,pos,deltahi);
135     vectorAbs3D(deltalo);
136     vectorAbs3D(deltahi);
137 
138     dlo = vectorMin3D(deltalo);
139     dhi = vectorMin3D(deltahi);
140 
141     if(dlo < dhi)
142         return dlo;
143     return dhi;
144 }
145 
146 /* ----------------------------------------------------------------------
147    return smallest extent ob subbox
148 ------------------------------------------------------------------------- */
149 
min_subbox_extent(double & min_extent,int & dim)150 inline void Domain::min_subbox_extent(double &min_extent,int &dim)
151 {
152     if(is_wedge)
153         error->one(FLERR,"missing implementation");
154 
155     double delta[3];
156     vectorSubtract3D(subhi,sublo,delta);
157 
158     min_extent = vectorMin3D(delta,dim);
159 }
160 
161 /* ----------------------------------------------------------------------
162    domain check
163 ------------------------------------------------------------------------- */
164 
is_periodic_ghost(int i)165 inline int Domain::is_periodic_ghost(int i)
166 {
167     if(i < atom->nlocal) return 0;
168 
169     if(is_wedge)
170         return is_periodic_ghost_wedge(i);
171 
172     int idim;
173     double *x = atom->x[i];
174     const double cutneighmax = neighbor->cutneighmax;
175 
176     for(idim = 0; idim < 3; idim++)
177          if ((x[idim] < (boxlo[idim]+cutneighmax) || x[idim] > (boxhi[idim]-cutneighmax)) && periodicity[idim])
178 
179             return 1;
180 
181     return 0;
182 }
183 
184 /* ----------------------------------------------------------------------
185    check if atom is the true representation of the particle on this subdomain
186    used when tallying stats across owned and ghost particles
187 ------------------------------------------------------------------------- */
188 
is_owned_or_first_ghost(int i)189 inline bool Domain::is_owned_or_first_ghost(int i)
190 {
191     if(!atom->tag_enable)
192         error->one(FLERR,"The current simulation setup requires atoms to have tags");
193     if(0 == atom->map_style)
194         error->one(FLERR,"The current simulation setup requires an 'atom_modify map' command to allocate an atom map");
195 
196     if(i == atom->map(atom->tag[i]))
197         return true;
198     return false;
199 }
200 
201 #endif
202