1 #include "PerPairQuantity.h"
2 #include "PerAtomQuantity.h"
3 #include "KernelFunction.h"
4 #include "FE_Mesh.h"
5 #include "Utility.h"
6 #include "Quadrature.h"
7
8 using ATC::HeartBeat;
9 using std::pair;
10 using std::map;
11
12 namespace ATC {
13
14 //==========================================================
PairMap(LammpsInterface * lammpsInterface,int groupbit)15 PairMap::PairMap(LammpsInterface * lammpsInterface, int groupbit ):
16 lammpsInterface_(lammpsInterface),
17 groupbit_(groupbit),
18 nPairs_(0), nBonds_(0)
19 {
20 };
~PairMap(void)21 PairMap::~PairMap(void)
22 {
23 };
24 //==========================================================
PairMapNeighbor(LammpsInterface * lammpsInterface,int groupbit)25 PairMapNeighbor::PairMapNeighbor(LammpsInterface * lammpsInterface, int groupbit):
26 PairMap(lammpsInterface,groupbit)
27 {
28 };
29
reset(void) const30 void PairMapNeighbor::reset(void) const
31 {
32 int inum = lammpsInterface_->neighbor_list_inum();
33 int *ilist = lammpsInterface_->neighbor_list_ilist();
34 int *numneigh = lammpsInterface_->neighbor_list_numneigh();
35 int **firstneigh = lammpsInterface_->neighbor_list_firstneigh();
36 const int * mask = lammpsInterface_->atom_mask();
37
38 pairMap_.clear();
39 int pairIndex = nBonds_;
40 std::pair< int,int > pair_ij;
41 for (int i = 0; i < inum; i++) {
42 int lammps_i = ilist[i];
43 if (mask[lammps_i] & groupbit_) {
44 for (int j = 0; j < numneigh[lammps_i]; j++) {
45 int lammps_j = firstneigh[lammps_i][j];
46 lammpsInterface_->neighbor_remap(lammps_j);
47 pair_ij.first = lammps_i; // alpha
48 pair_ij.second = lammps_j; // beta
49 pairMap_[pair_ij] = pairIndex;
50 pairIndex++;
51 }
52 }
53 }
54 nPairs_ = pairIndex;
55 needReset_ = false;
56 }
57
58 //==========================================================
PairMapBond(LammpsInterface * lammpsInterface,int groupbit)59 PairMapBond::PairMapBond(LammpsInterface * lammpsInterface, int groupbit):
60 PairMap(lammpsInterface,groupbit)
61 {
62 };
63
64 //==========================================================
PairMapBoth(LammpsInterface * lammpsInterface,int groupbit)65 PairMapBoth::PairMapBoth(LammpsInterface * lammpsInterface, int groupbit):
66 PairMapNeighbor(lammpsInterface,groupbit)
67 {
68 };
69 //==========================================================
DensePerPairMatrix(LammpsInterface * lammpsInterface,const PairMap & pairMap,int nCols)70 DensePerPairMatrix::DensePerPairMatrix(LammpsInterface * lammpsInterface,
71 const PairMap & pairMap,
72 int nCols):
73 lammpsInterface_(lammpsInterface),
74 pairMap_(pairMap),
75 nCols_(nCols)
76 {
77 };
78
79 //==========================================================
PairVirial(LammpsInterface * lammpsInterface,const PairMap & pairMap,int nCols)80 PairVirial::PairVirial(LammpsInterface * lammpsInterface,
81 const PairMap & pairMap, int nCols):
82 DensePerPairMatrix(lammpsInterface,pairMap,nCols)
83 {
84 };
85 //==========================================================
PairVirialEulerian(LammpsInterface * lammpsInterface,const PairMap & pairMap)86 PairVirialEulerian::PairVirialEulerian(LammpsInterface * lammpsInterface,
87 const PairMap & pairMap):
88 PairVirial(lammpsInterface,pairMap,6)
89 {
90 };
91
92
reset(void) const93 void PairVirialEulerian::reset(void) const
94 {
95 int nPairs = pairMap_.size();
96 quantity_.reset(nPairs,nCols_);
97
98 double **xatom = lammpsInterface_->xatom();
99 for (ATOM_PAIR apair = pairMap_.start();
100 ! pairMap_.finished(); apair=pairMap_++){
101 int lammps_a = (apair.first).first ;
102 int lammps_b = (apair.first).second;
103 int pairIndex = apair.second;
104 double * xa = xatom[lammps_a];
105 double * xb = xatom[lammps_b];
106 double delx = xa[0] - xb[0];
107 double dely = xa[1] - xb[1];
108 double delz = xa[2] - xb[2];
109 double rsq = delx*delx + dely*dely + delz*delz;
110 double fforce = 0;
111 lammpsInterface_->pair_force(apair,rsq,fforce);
112 quantity_(pairIndex,0)=-delx*delx*fforce;
113 quantity_(pairIndex,1)=-dely*dely*fforce;
114 quantity_(pairIndex,2)=-delz*delz*fforce;
115 quantity_(pairIndex,3)=-delx*dely*fforce;
116 quantity_(pairIndex,4)=-delx*delz*fforce;
117 quantity_(pairIndex,5)=-dely*delz*fforce;
118 }
119 }
120 //==========================================================
PairVirialLagrangian(LammpsInterface * lammpsInterface,const PairMap & pairMap,double ** xRef)121 PairVirialLagrangian::PairVirialLagrangian(LammpsInterface * lammpsInterface,
122 const PairMap & pairMap,
123 double ** xRef):
124 // const PerAtomQuantity<double> * xRef):
125 PairVirial(lammpsInterface,pairMap,9),
126 xRef_(xRef)
127 {
128
129 };
130
131
reset(void) const132 void PairVirialLagrangian::reset(void) const
133 {
134 int nPairs = pairMap_.size();
135 quantity_.reset(nPairs,nCols_);
136
137 double **xatom = lammpsInterface_->xatom();
138
139 double ** xref = xRef_;
140
141 for (ATOM_PAIR apair = pairMap_.start();
142 ! pairMap_.finished(); apair=pairMap_++){
143 int lammps_a = (apair.first).first ;
144 int lammps_b = (apair.first).second;
145 int pairIndex = apair.second;
146 double * xa = xatom[lammps_a];
147 double * xb = xatom[lammps_b];
148 double delx = xa[0] - xb[0];
149 double dely = xa[1] - xb[1];
150 double delz = xa[2] - xb[2];
151 double * Xa = xref[lammps_a];
152 double * Xb = xref[lammps_b];
153 double delX = Xa[0] - Xb[0];
154 double delY = Xa[1] - Xb[1];
155 double delZ = Xa[2] - Xb[2];
156 double rsq = delx*delx + dely*dely + delz*delz;
157 double fforce = 0;
158 lammpsInterface_->pair_force(apair,rsq,fforce);
159 quantity_(pairIndex,0)=-delx*fforce*delX;
160 quantity_(pairIndex,1)=-delx*fforce*delY;
161 quantity_(pairIndex,2)=-delx*fforce*delZ;
162 quantity_(pairIndex,3)=-dely*fforce*delX;
163 quantity_(pairIndex,4)=-dely*fforce*delY;
164 quantity_(pairIndex,5)=-dely*fforce*delZ;
165 quantity_(pairIndex,6)=-delz*fforce*delX;
166 quantity_(pairIndex,7)=-delz*fforce*delY;
167 quantity_(pairIndex,8)=-delz*fforce*delZ;
168 }
169 }
170 //==========================================================
PairPotentialHeatFlux(LammpsInterface * lammpsInterface,const PairMap & pairMap)171 PairPotentialHeatFlux::PairPotentialHeatFlux(LammpsInterface * lammpsInterface,
172 const PairMap & pairMap):
173 DensePerPairMatrix(lammpsInterface,pairMap,3)
174 {
175 };
176 //==========================================================
PairPotentialHeatFluxEulerian(LammpsInterface * lammpsInterface,const PairMap & pairMap)177 PairPotentialHeatFluxEulerian::PairPotentialHeatFluxEulerian(LammpsInterface * lammpsInterface,
178 const PairMap & pairMap):
179 PairPotentialHeatFlux(lammpsInterface,pairMap)
180 {
181
182 };
183
reset(void) const184 void PairPotentialHeatFluxEulerian::reset(void) const
185 {
186 int nPairs = pairMap_.size();
187 quantity_.reset(nPairs,nCols_);
188
189 double **xatom = lammpsInterface_->xatom();
190 double **vatom = lammpsInterface_->vatom();
191 for (ATOM_PAIR apair = pairMap_.start();
192 ! pairMap_.finished(); apair=pairMap_++){
193 int lammps_a = (apair.first).first ;
194 int lammps_b = (apair.first).second;
195 int pairIndex = apair.second;
196 double * xa = xatom[lammps_a];
197 double * xb = xatom[lammps_b];
198 double delx = xa[0] - xb[0];
199 double dely = xa[1] - xb[1];
200 double delz = xa[2] - xb[2];
201 double rsq = delx*delx + dely*dely + delz*delz;
202 double fforce = 0;
203 lammpsInterface_->pair_force(apair,rsq,fforce);
204 double* v = vatom[lammps_a];
205 fforce *=delx*v[0] + dely*v[1] + delz*v[2];
206 quantity_(pairIndex,0)=fforce*delx;
207 quantity_(pairIndex,1)=fforce*dely;
208 quantity_(pairIndex,2)=fforce*delz;
209 }
210 }
211 //==========================================================
PairPotentialHeatFluxLagrangian(LammpsInterface * lammpsInterface,const PairMap & pairMap,double ** xRef)212 PairPotentialHeatFluxLagrangian::PairPotentialHeatFluxLagrangian(LammpsInterface * lammpsInterface,
213 const PairMap & pairMap, double ** xRef):
214 PairPotentialHeatFlux(lammpsInterface,pairMap),
215 xRef_(xRef)
216 {
217
218 };
219
reset(void) const220 void PairPotentialHeatFluxLagrangian::reset(void) const
221 {
222 int nPairs = pairMap_.size();
223 quantity_.reset(nPairs,nCols_);
224
225 double **xatom = lammpsInterface_->xatom();
226 double **vatom = lammpsInterface_->vatom();
227
228 for (ATOM_PAIR apair = pairMap_.start();
229 ! pairMap_.finished(); apair=pairMap_++){
230 int lammps_a = (apair.first).first ;
231 int lammps_b = (apair.first).second;
232 int pairIndex = apair.second;
233 double * xa = xatom[lammps_a];
234 double * xb = xatom[lammps_b];
235 double delx = xa[0] - xb[0];
236 double dely = xa[1] - xb[1];
237 double delz = xa[2] - xb[2];
238 double * Xa = xRef_[lammps_a];
239 double * Xb = xRef_[lammps_b];
240 double delX = Xa[0] - Xb[0];
241 double delY = Xa[1] - Xb[1];
242 double delZ = Xa[2] - Xb[2];
243 double rsq = delx*delx + dely*dely + delz*delz;
244 double fforce = 0;
245 lammpsInterface_->pair_force(apair,rsq,fforce);
246 double* v = vatom[lammps_a];
247 fforce *=delx*v[0] + dely*v[1] + delz*v[2];
248 quantity_(pairIndex,0)=fforce*delX;
249 quantity_(pairIndex,1)=fforce*delY;
250 quantity_(pairIndex,2)=fforce*delZ;
251 }
252 }
253 //==========================================================
SparsePerPairMatrix(LammpsInterface * lammpsInterface,const PairMap & pairMap)254 SparsePerPairMatrix::SparsePerPairMatrix(LammpsInterface * lammpsInterface,
255 const PairMap & pairMap):
256 lammpsInterface_(lammpsInterface),
257 pairMap_(pairMap)
258 {
259 };
260 //==========================================================
BondMatrix(LammpsInterface * lammpsInterface,const PairMap & pairMap,double ** x,const FE_Mesh * feMesh)261 BondMatrix::BondMatrix(LammpsInterface * lammpsInterface,
262 const PairMap & pairMap, double ** x, const FE_Mesh * feMesh):
263 SparsePerPairMatrix(lammpsInterface,pairMap), x_(x), feMesh_(feMesh)
264 {
265 };
266 //==========================================================
BondMatrixKernel(LammpsInterface * lammpsInterface,const PairMap & pairMap,double ** x,const FE_Mesh * feMesh,const KernelFunction * kernelFunction)267 BondMatrixKernel::BondMatrixKernel(LammpsInterface * lammpsInterface,
268 const PairMap & pairMap,
269 double ** x,
270 const FE_Mesh * feMesh,
271 const KernelFunction * kernelFunction):
272 BondMatrix(lammpsInterface,pairMap,x,feMesh),
273 kernelFunction_(kernelFunction)
274 {
275 };
reset(void) const276 void BondMatrixKernel::reset(void) const
277 {
278 int nPairs = pairMap_.size(); // needs to come after quantity for reset
279 int nNodes = feMesh_->num_nodes_unique();
280 quantity_.reset(nNodes,nPairs);
281 double lam1,lam2;
282 std::pair< int,int > pair_jk;
283 int heartbeatFreq = (nNodes <= 10 ? 1 : (int) nNodes / 10);
284 HeartBeat beat("computing bond matrix ",heartbeatFreq);
285 beat.start();
286 DENS_VEC xa(3),xI(3),xaI(3),xb(3),xbI(3),xba(3);
287 double invVol = kernelFunction_->inv_vol();
288 for (int I = 0; I < nNodes; I++) {
289 beat.next();
290 xI = feMesh_->nodal_coordinates(I);
291 if (!kernelFunction_->node_contributes(xI)) { continue; }
292 for (ATOM_PAIR apair = pairMap_.start();
293 ! pairMap_.finished(); apair=pairMap_++){
294 int lammps_a = (apair.first).first ;
295 int lammps_b = (apair.first).second;
296 xa.copy(x_[lammps_a],3);
297 xaI = xa - xI;
298 lammpsInterface_->periodicity_correction(xaI.ptr());
299 xb.copy(x_[lammps_b],3);
300 xba = xb - xa;
301 xbI = xba + xaI;
302 kernelFunction_->bond_intercepts(xaI,xbI,lam1,lam2);
303 if (lam1 < lam2) {
304 double bondValue = invVol*(kernelFunction_->bond(xaI,xbI,lam1,lam2));
305 int pairIndex = apair.second;
306 quantity_.add(I,pairIndex,bondValue);
307 } // if lam1 < lam2
308 } // pair map
309 } // end nodes loop
310 quantity_.compress();
311 beat.finish();
312 }
313 //==========================================================
BondMatrixPartitionOfUnity(LammpsInterface * lammpsInterface,const PairMap & pairMap,double ** x,const FE_Mesh * feMesh,const DIAG_MAN * invVols)314 BondMatrixPartitionOfUnity::BondMatrixPartitionOfUnity(LammpsInterface * lammpsInterface,
315 const PairMap & pairMap, double ** x, const FE_Mesh * feMesh,
316 const DIAG_MAN * invVols):
317 BondMatrix(lammpsInterface,pairMap,x,feMesh),
318 invVols_(invVols)
319 {
320 ATC::Quadrature::instance()->set_line_quadrature(lineNgauss_,lineXg_,lineWg_);
321 double lam1 = 0.0, lam2 = 1.0;
322 double del_lambda = 0.5*(lam2 - lam1);
323 double avg_lambda = 0.5*(lam2 + lam1);
324 for (int i = 0; i < lineNgauss_; i++) {
325 double lambda = del_lambda*lineXg_[i] +avg_lambda;
326 lineXg_[i] = lambda;
327 lineWg_[i] *= 0.5;
328 }
329 };
reset(void) const330 void BondMatrixPartitionOfUnity::reset(void) const
331 {
332 int nNodes = feMesh_->num_nodes_unique();
333 int nPairs = pairMap_.size();
334 quantity_.reset(nNodes,nPairs);
335 int nodes_per_element = feMesh_->num_nodes_per_element();
336 Array<int> node_list(nodes_per_element);
337 DENS_VEC shp(nodes_per_element);
338 std::pair< int,int > pair_jk;
339 int heartbeatFreq = (int) nPairs / 10;
340 HeartBeat beat("computing bond matrix ",heartbeatFreq);
341 beat.start();
342 DENS_VEC xa(3),xb(3),xab(3),xlambda(3);
343 for (ATOM_PAIR apair = pairMap_.start();
344 ! pairMap_.finished(); apair=pairMap_++){
345 beat.next();
346 int lammps_a = (apair.first).first ;
347 int lammps_b = (apair.first).second;
348 int pairIndex = apair.second;
349 xa.copy(x_[lammps_a],3);
350 xb.copy(x_[lammps_b],3);
351 xab = xa - xb;
352 for (int i = 0; i < lineNgauss_; i++) {
353 double lambda = lineXg_[i];
354 xlambda = lambda*xab + xb;
355 lammpsInterface_->periodicity_correction(xlambda.ptr());
356 feMesh_->shape_functions(xlambda,shp,node_list);
357 // accumulate to nodes whose support overlaps the integration point
358 for (int I = 0; I < nodes_per_element; I++) {
359 int Inode = node_list(I);
360 double inv_vol = (invVols_->quantity())(Inode,Inode);
361 double val = inv_vol*shp(I)*lineWg_[i];
362 quantity_.add(Inode,pairIndex,val);
363 }
364 }
365 }
366 quantity_.compress();
367 beat.finish();
368 }
369 //==========================================================
370 }
371