1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#include "region_block_kokkos.h"
#include "atom_kokkos.h"
#include "atom_masks.h"
using namespace LAMMPS_NS;
#define BIG 1.0e20
/* ---------------------------------------------------------------------- */
template<class DeviceType>
RegBlockKokkos<DeviceType>::RegBlockKokkos(LAMMPS *lmp, int narg, char **arg) : RegBlock(lmp, narg, arg)
{
atomKK = (AtomKokkos*) atom;
}
/* ---------------------------------------------------------------------- */
template<class DeviceType>
RegBlockKokkos<DeviceType>::~RegBlockKokkos()
{
}
/* ----------------------------------------------------------------------
inside = 1 if x,y,z is inside or on surface
inside = 0 if x,y,z is outside and not on surface
------------------------------------------------------------------------- */
template<class DeviceType>
KOKKOS_INLINE_FUNCTION
int RegBlockKokkos<DeviceType>::k_inside(double x, double y, double z) const
{
if (x >= xlo && x <= xhi && y >= ylo && y <= yhi && z >= zlo && z <= zhi)
return 1;
return 0;
}
template<class DeviceType>
void RegBlockKokkos<DeviceType>::match_all_kokkos(int groupbit_in, DAT::tdual_int_1d k_match_in)
{
groupbit = groupbit_in;
d_match = k_match_in.template view<DeviceType>();
atomKK->sync(Device, X_MASK | MASK_MASK);
x = atomKK->k_x.view<DeviceType>();
mask = atomKK->k_mask.view<DeviceType>();
int nlocal = atom->nlocal;
copymode = 1;
Kokkos::parallel_for(Kokkos::RangePolicy<DeviceType, TagRegBlockMatchAll>(0,nlocal),*this);
copymode = 0;
k_match_in.template modify<DeviceType>();
}
template<class DeviceType>
KOKKOS_INLINE_FUNCTION
void RegBlockKokkos<DeviceType>::operator()(TagRegBlockMatchAll, const int &i) const {
if (mask[i] & groupbit) {
double x_tmp = x(i,0);
double y_tmp = x(i,1);
double z_tmp = x(i,2);
d_match[i] = match(x_tmp,y_tmp,z_tmp);
}
}
/* ----------------------------------------------------------------------
determine if point x,y,z is a match to region volume
XOR computes 0 if 2 args are the same, 1 if different
note that k_inside() returns 1 for points on surface of region
thus point on surface of exterior region will not match
if region has variable shape, invoke shape_update() once per timestep
if region is dynamic, apply inverse transform to x,y,z
unmove first, then unrotate, so don't have to change rotation point
caller is responsible for wrapping this call with
modify->clearstep_compute() and modify->addstep_compute() if needed
------------------------------------------------------------------------- */
template<class DeviceType>
KOKKOS_INLINE_FUNCTION
int RegBlockKokkos<DeviceType>::match(double x, double y, double z) const
{
if (dynamic) inverse_transform(x,y,z);
return !(k_inside(x,y,z) ^ interior);
}
/* ----------------------------------------------------------------------
transform a point x,y,z in moved space back to region space
undisplace first, then unrotate (around original P)
------------------------------------------------------------------------- */
template<class DeviceType>
KOKKOS_INLINE_FUNCTION
void RegBlockKokkos<DeviceType>::inverse_transform(double &x, double &y, double &z) const
{
if (moveflag) {
x -= dx;
y -= dy;
z -= dz;
}
if (rotateflag) rotate(x,y,z,-theta);
}
/* ----------------------------------------------------------------------
rotate x,y,z by angle via right-hand rule around point and runit normal
sign of angle determines whether rotating forward/backward in time
return updated x,y,z
R = vector axis of rotation
P = point = point to rotate around
R0 = runit = unit vector for R
X0 = x,y,z = initial coord of atom
D = X0 - P = vector from P to X0
C = (D dot R0) R0 = projection of D onto R, i.e. Dparallel
A = D - C = vector from R line to X0, i.e. Dperp
B = R0 cross A = vector perp to A in plane of rotation, same len as A
A,B define plane of circular rotation around R line
new x,y,z = P + C + A cos(angle) + B sin(angle)
------------------------------------------------------------------------- */
template<class DeviceType>
KOKKOS_INLINE_FUNCTION
void RegBlockKokkos<DeviceType>::rotate(double &x, double &y, double &z, double angle) const
{
double a[3],b[3],c[3],d[3],disp[3];
double sine = sin(angle);
double cosine = cos(angle);
d[0] = x - point[0];
d[1] = y - point[1];
d[2] = z - point[2];
double x0dotr = d[0]*runit[0] + d[1]*runit[1] + d[2]*runit[2];
c[0] = x0dotr * runit[0];
c[1] = x0dotr * runit[1];
c[2] = x0dotr * runit[2];
a[0] = d[0] - c[0];
a[1] = d[1] - c[1];
a[2] = d[2] - c[2];
b[0] = runit[1]*a[2] - runit[2]*a[1];
b[1] = runit[2]*a[0] - runit[0]*a[2];
b[2] = runit[0]*a[1] - runit[1]*a[0];
disp[0] = a[0]*cosine + b[0]*sine;
disp[1] = a[1]*cosine + b[1]*sine;
disp[2] = a[2]*cosine + b[2]*sine;
x = point[0] + c[0] + disp[0];
y = point[1] + c[1] + disp[1];
z = point[2] + c[2] + disp[2];
}
namespace LAMMPS_NS {
template class RegBlockKokkos<LMPDeviceType>;
#ifdef KOKKOS_ENABLE_CUDA
template class RegBlockKokkos<LMPHostType>;
#endif
}