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