MDStressLab++
MethodLdad.cpp
Go to the documentation of this file.
1 /*
2  * MethodLdad.cpp
3  *
4  * Created on: Nov 5, 2019
5  * Author: Nikhil
6  */
7 
8 #include <math.h>
9 #include <iostream>
10 #include <limits>
11 #include "typedef.h"
12 
13 int PointLineRelationship(const double& p);
14 
15 template<typename T>
16 MethodLdad<T>::MethodLdad(const Matrix3d& ldadVectors): ldadVectors(ldadVectors)
17 {
18  // initialize the normalizer here using oneDFunction.integrate(-1,1)
19  normalizer = 1.0/(8.0*fabs(ldadVectors.determinant())) ;
20 
21  // initialize the averagingDomainSize
22  double p1, p2, p3;
23  p1 = fabs(ldadVectors(0,0)) + fabs(ldadVectors(0,1)) + fabs(ldadVectors(0,2));
24  p2 = fabs(ldadVectors(1,0)) + fabs(ldadVectors(1,1)) + fabs(ldadVectors(1,2));
25  p3 = fabs(ldadVectors(2,0)) + fabs(ldadVectors(2,1)) + fabs(ldadVectors(2,2));
26 
27 
28  this->averagingDomainSize = sqrt(p1 * p1 + p2 * p2 + p3 * p3);
29 
30  // initialize the inverseLdadVectors
31  inverseLdadVectors = ldadVectors.inverse();
32 }
33 
34 // template<typename T>
35 // MethodLdad<T>::MethodLdad(const MethodLdad<T>& _MethodLdad)
36 // {
37 // *this= _MethodLdad;
38 // }
39 
40 template<typename T>
41 MethodLdad<T>::~MethodLdad() {
42  // TODO Auto-generated destructor stub
43 }
44 
45 template<typename T>
46 double MethodLdad<T>::operator()(const Vector3d& vec) const
47 {
48  Vector3d vec_pull;
49  vec_pull = inverseLdadVectors * vec.transpose();
50  // oneDFunction -1 1
51  return oneDFunction(vec_pull(0))*oneDFunction(vec_pull(1))*oneDFunction(vec_pull(2));
52 }
53 
54 template<typename T>
55 double MethodLdad<T>::bondFunction(const Vector3d& vec1, const Vector3d& vec2) const
56 {
57  /* copy from fortran
58  The algorithm is:
59  selecting two points out of eight points.
60  The eight points composed of two points of vec1 and vec2 themselves,
61  and the intersection of x = +-1, y = +-1, z = +-1. */
62 
63  Vector3d vec1_pull, vec2_pull, vec12_pull;
64  Vector3d vec1_pull_seg, vec2_pull_seg;
65  Vector3d vec1_push_seg, vec2_push_seg, vec12_push_seg;
66 
67  Vector3i degenerate(0,0,0);
68  VectorXi selected(8);
69  int selected_count = 0;
70  Matrix3i relative_position;
71 
72  ArrayXXd r_pull_intersect(8,3);
73  double total_length = 0.0;
74  double distance = (vec2 - vec1).norm();
75 
76  for (int i = 0; i < 8;i++)
77  {
78  for (int j = 0; j < 3; j++)
79  {
80  r_pull_intersect(i,j) = 0.0;
81  }
82  selected(i) = 0;
83  }
84 
85  vec1_pull = inverseLdadVectors * vec1.transpose();
86  vec2_pull = inverseLdadVectors * vec2.transpose();
87  vec12_pull = vec2_pull - vec1_pull;
88  // This is a necessary step to initialization
89  vec1_pull_seg = vec1_pull;
90  vec2_pull_seg = vec2_pull;
91 
92  for (int i = 0; i <= 2; i++)
93  {
94  if (fabs(vec12_pull(i)) < epsilon)
95  {
96  degenerate(i) = 1;
97  }
98  else
99  {
100  degenerate(i) = 0;
101  }
102  }
103 
104  relative_position(0,0) = PointLineRelationship(vec1_pull(0));
105  relative_position(1,0) = PointLineRelationship(vec1_pull(1));
106  relative_position(2,0) = PointLineRelationship(vec1_pull(2));
107  relative_position(0,1) = PointLineRelationship(vec2_pull(0));
108  relative_position(1,1) = PointLineRelationship(vec2_pull(1));
109  relative_position(2,1) = PointLineRelationship(vec2_pull(2));
110  relative_position(0,2) = 0;
111  relative_position(1,2) = 0;
112  relative_position(2,2) = 0;
113 
114  // degenerate to 0D. It is impossible
115  if (degenerate(0) && degenerate(1) && degenerate(2))
116  {
117  MY_ERROR("Degeneration to 0D. This indicates you have overlapping atoms in your system!");
118  }
119  // degenerate to 1D x
120  else if (!degenerate(0) && degenerate(1) && degenerate(2))
121  {
122  relative_position(1,2) = PointLineRelationship((vec1_pull(1) + vec2_pull(1)) / 2.0);
123  relative_position(2,2) = PointLineRelationship((vec1_pull(2) + vec2_pull(2)) / 2.0);
124  if (relative_position(1,2) == 1 || relative_position(1,2) == 5 || \
125  relative_position(2,2) == 1 || relative_position(2,2) == 5)
126  {
127  return 0.0;
128  }
129  // determine x length
130  if (relative_position(0,0) == 2 || relative_position(0,0) == 3 || relative_position(0,0) == 4)
131  {
132  vec1_pull_seg(0) = vec1_pull(0);
133  if (relative_position(0,1) == 2 || relative_position(0,1) == 3 || relative_position(0,1) == 4)
134  {
135  vec2_pull_seg(0) = vec2_pull(0);
136  }
137  else if (relative_position(0,1) == 1)
138  {
139  vec2_pull_seg(0) = -1.0;
140  }
141  else if (relative_position(0,1) == 5)
142  {
143  vec2_pull_seg(0) = 1.0;
144  }
145  }
146  else if (relative_position(0,0) == 1)
147  {
148  if (relative_position(0,1) == 2 || relative_position(0,1) == 3 || relative_position(0,1) == 4)
149  {
150  vec1_pull_seg(0) = -1.0;
151  vec2_pull_seg(0) = vec2_pull(0);
152  }
153  else if (relative_position(0,1) == 1)
154  {
155  return 0.0;
156  }
157  else if (relative_position(0,1) == 5)
158  {
159  vec1_pull_seg(0) = -1.0;
160  vec2_pull_seg(0) = 1.0;
161  }
162  }
163  else if (relative_position(0,0) == 5)
164  {
165  if (relative_position(0,1) == 2 || relative_position(0,1) == 3 || relative_position(0,1) == 4)
166  {
167  vec1_pull_seg(0) = 1.0;
168  vec2_pull_seg(0) = vec2_pull(0);
169  }
170  else if (relative_position(0,1) == 1)
171  {
172  vec1_pull_seg(0) = 1.0;
173  vec2_pull_seg(0) = -1.0;
174  }
175  else if (relative_position(0,1) == 5)
176  {
177  return 0.0;
178  }
179  }
180  }
181  // degenerate to 1D y
182  else if (degenerate(0) && !degenerate(1) && degenerate(2))
183  {
184  relative_position(0,2) = PointLineRelationship((vec1_pull(0) + vec2_pull(0)) / 2.0);
185  relative_position(2,2) = PointLineRelationship((vec1_pull(2) + vec2_pull(2)) / 2.0);
186  if (relative_position(0,2) == 1 || relative_position(0,2) == 5 || \
187  relative_position(2,2) == 1 || relative_position(2,2) == 5)
188  {
189  return 0.0;
190  }
191  // determine y length
192  if (relative_position(1,0) == 2 || relative_position(1,0) == 3 || relative_position(1,0) == 4)
193  {
194  vec1_pull_seg(1) = vec1_pull(1);
195  if (relative_position(1,1) == 2 || relative_position(1,1) == 3 || relative_position(1,1) == 4)
196  {
197  vec2_pull_seg(1) = vec2_pull(1);
198  }
199  else if (relative_position(1,1) == 1)
200  {
201  vec2_pull_seg(1) = -1.0;
202  }
203  else if (relative_position(1,1) == 5)
204  {
205  vec2_pull_seg(1) = 1.0;
206  }
207  }
208  else if (relative_position(1,0) == 1)
209  {
210  if (relative_position(1,1) == 2 || relative_position(1,1) == 3 || relative_position(1,1) == 4)
211  {
212  vec1_pull_seg(1) = -1.0;
213  vec2_pull_seg(1) = vec2_pull(1);
214  }
215  else if (relative_position(1,1) == 1)
216  {
217  return 0.0;
218  }
219  else if (relative_position(1,1) == 5)
220  {
221  vec1_pull_seg(1) = -1.0;
222  vec2_pull_seg(1) = 1.0;
223  }
224  }
225  else if (relative_position(1,0) == 5)
226  {
227  if (relative_position(1,1) == 2 || relative_position(1,1) == 3 || relative_position(1,1) == 4)
228  {
229  vec1_pull_seg(1) = 1.0;
230  vec2_pull_seg(1) = vec2_pull(1);
231  }
232  else if (relative_position(1,1) == 1)
233  {
234  vec1_pull_seg(1) = 1.0;
235  vec2_pull_seg(1) = -1.0;
236  }
237  else if (relative_position(1,1) == 5)
238  {
239  return 0.0;
240  }
241  }
242  }
243  // degenerate to 1D z
244  else if (degenerate(0) && degenerate(1) && !degenerate(2))
245  {
246  relative_position(0,2) = PointLineRelationship((vec1_pull(0) + vec2_pull(0)) / 2.0);
247  relative_position(1,2) = PointLineRelationship((vec1_pull(1) + vec2_pull(1)) / 2.0);
248  if (relative_position(0,2) == 1 || relative_position(0,2) == 5 || \
249  relative_position(1,2) == 1 || relative_position(1,2) == 5)
250  {
251  return 0.0;
252  }
253  // determine z length
254  if (relative_position(2,0) == 2 || relative_position(2,0) == 3 || relative_position(2,0) == 4)
255  {
256  vec1_pull_seg(2) = vec1_pull(2);
257  if (relative_position(2,1) == 2 || relative_position(2,1) == 3 || relative_position(2,1) == 4)
258  {
259  vec2_pull_seg(2) = vec2_pull(2);
260  }
261  else if (relative_position(2,1) == 1)
262  {
263  vec2_pull_seg(2) = -1.0;
264  }
265  else if (relative_position(2,1) == 5)
266  {
267  vec2_pull_seg(2) = 1.0;
268  }
269  }
270  else if (relative_position(2,0) == 1)
271  {
272  if (relative_position(2,1) == 2 || relative_position(2,1) == 3 || relative_position(2,1) == 4)
273  {
274  vec1_pull_seg(2) = -1.0;
275  vec2_pull_seg(2) = vec2_pull(2);
276  }
277  else if (relative_position(2,1) == 1)
278  {
279  return 0.0;
280  }
281  else if (relative_position(2,1) == 5)
282  {
283  vec1_pull_seg(2) = -1.0;
284  vec2_pull_seg(2) = 1.0;
285  }
286  }
287  else if (relative_position(2,0) == 5)
288  {
289  if (relative_position(2,1) == 2 || relative_position(2,1) == 3 || relative_position(2,1) == 4)
290  {
291  vec1_pull_seg(2) = 1.0;
292  vec2_pull_seg(2) = vec2_pull(2);
293  }
294  else if (relative_position(2,1) == 1)
295  {
296  vec1_pull_seg(2) = 1.0;
297  vec2_pull_seg(2) = -1.0;
298  }
299  else if (relative_position(2,1) == 5)
300  {
301  return 0.0;
302  }
303  }
304  }
305  // degenerate to 2D xy
306  else if (!degenerate(0) && !degenerate(1) && degenerate(2))
307  {
308  relative_position(2,2) = PointLineRelationship((vec1_pull(2) + vec2_pull(2)) / 2.0);
309  if (relative_position(2,2) == 1 || relative_position(2,2) == 5)
310  {
311  return 0.0;
312  }
313 
314  for (int i = 0; i < 8; i++)
315  {
316  selected(i) = 0;
317  }
318 
319  // x = -1.0
320  r_pull_intersect(0,0) = -1.0;
321  r_pull_intersect(0,1) = vec1_pull(1) + \
322  (vec2_pull(1) - vec1_pull(1)) / (vec2_pull(0) - vec1_pull(0)) * (-1.0 - vec1_pull(0));
323  // x = 1.0
324  r_pull_intersect(1,0) = 1.0;
325  r_pull_intersect(1,1) = vec1_pull(1) + \
326  (vec2_pull(1) - vec1_pull(1)) / (vec2_pull(0) - vec1_pull(0)) * (1.0 - vec1_pull(0));
327  // y = -1.0
328  r_pull_intersect(2,0) = vec1_pull(0) + \
329  (vec2_pull(0) - vec1_pull(0)) / (vec2_pull(1) - vec1_pull(1)) * (-1.0 - vec1_pull(1));
330  r_pull_intersect(2,1) = -1.0;
331  // y = 1.0
332  r_pull_intersect(3,0) = vec1_pull(0) + \
333  (vec2_pull(0) - vec1_pull(0)) / (vec2_pull(1) - vec1_pull(1)) * (1.0 - vec1_pull(1));
334  r_pull_intersect(3,1) = 1.0;
335 
336  // vec1 itself
337  r_pull_intersect(6,0) = vec1_pull(0);
338  r_pull_intersect(6,1) = vec1_pull(1);
339  r_pull_intersect(6,2) = vec1_pull(2);
340 
341  // vec2 itself
342  r_pull_intersect(7,0) = vec2_pull(0);
343  r_pull_intersect(7,1) = vec2_pull(1);
344  r_pull_intersect(7,2) = vec2_pull(2);
345 
346  // if the point is inside, then it is the point selected.
347  if ((relative_position(0,0) == 2 || relative_position(0,0) == 3 || relative_position(0,0) == 4) && \
348  (relative_position(1,0) == 2 || relative_position(1,0) == 3 || relative_position(1,0) == 4))
349  {
350  selected(6) = 1;
351  }
352 
353  if ((relative_position(0,1) == 2 || relative_position(0,1) == 3 || relative_position(0,1) == 4) && \
354  (relative_position(1,1) == 2 || relative_position(1,1) == 3 || relative_position(1,1) == 4))
355  {
356  selected(7) = 1;
357  }
358 
359  // must be inside cubes to be selected
360  for (int i = 0; i <= 1; i++)
361  {
362  if (((r_pull_intersect(i,0) > vec1_pull(0) && r_pull_intersect(i,0) < vec2_pull(0)) || \
363  (r_pull_intersect(i,0) > vec2_pull(0) && r_pull_intersect(i,0) < vec1_pull(0))) && \
364  ((r_pull_intersect(i,1) > vec1_pull(1) && r_pull_intersect(i,1) < vec2_pull(1)) || \
365  (r_pull_intersect(i,1) > vec2_pull(1) && r_pull_intersect(i,1) < vec1_pull(1))) && \
366  (r_pull_intersect(i,1) >= -1.0 - epsilon && r_pull_intersect(i,1) <= 1.0 + epsilon))
367  {
368  selected(i) = 1;
369  }
370  }
371  for (int i = 2; i <= 3; i++)
372  {
373  if (((r_pull_intersect(i,0) > vec1_pull(0) && r_pull_intersect(i,0) < vec2_pull(0)) || \
374  (r_pull_intersect(i,0) > vec2_pull(0) && r_pull_intersect(i,0) < vec1_pull(0))) && \
375  ((r_pull_intersect(i,1) > vec1_pull(1) && r_pull_intersect(i,1) < vec2_pull(1)) || \
376  (r_pull_intersect(i,1) > vec2_pull(1) && r_pull_intersect(i,1) < vec1_pull(1))) && \
377  (r_pull_intersect(i,0) >= -1.0 - epsilon && r_pull_intersect(i,0) <= 1.0 + epsilon))
378  {
379  selected(i) = 1;
380  }
381  }
382 
383  selected_count = 0;
384  for (int i = 0; i < 8; i++)
385  {
386  if (selected(i) && selected_count == 0)
387  {
388  vec1_pull_seg(0) = r_pull_intersect(i,0);
389  vec1_pull_seg(1) = r_pull_intersect(i,1);
390  selected_count = 1;
391  continue;
392  }
393  if (selected(i) && selected_count == 1)
394  {
395  if (fabs(vec1_pull_seg(0) - r_pull_intersect(i,0)) < epsilon && \
396  fabs(vec1_pull_seg(1) - r_pull_intersect(i,1)) < epsilon)
397  {
398  continue;
399  }
400  else
401  {
402  vec2_pull_seg(0) = r_pull_intersect(i,0);
403  vec2_pull_seg(1) = r_pull_intersect(i,1);
404  selected_count = 2;
405  continue;
406  }
407  }
408  if (selected(i) && selected_count == 2)
409  {
410  if ((fabs(vec1_pull_seg(0) - r_pull_intersect(i,0)) < epsilon && \
411  fabs(vec1_pull_seg(1) - r_pull_intersect(i,1)) < epsilon) || \
412  (fabs(vec2_pull_seg(0) - r_pull_intersect(i,0)) < epsilon && \
413  fabs(vec2_pull_seg(1) - r_pull_intersect(i,1)) < epsilon))
414  {
415  continue;
416  }
417  else
418  {
419  MY_ERROR("In LDAD bond_function xy, More than 2 different points \
420  with different positions are being selected. \
421  This means a line intersecting with a cube has 3 points. \
422  Something wrong with the algorithm. Need to print and debug!");
423  }
424  }
425  }
426  // No intersection .or. one point on the parallelepiped but another point is outside
427  if (selected_count == 0 || selected_count == 1)
428  {
429  return 0.0;
430  }
431  }
432  // degenerate to 2D xz
433  else if (!degenerate(0) && degenerate(1) && !degenerate(2))
434  {
435  relative_position(1,2) = PointLineRelationship((vec1_pull(1) + vec2_pull(1)) / 2.0);
436  if (relative_position(1,2) == 1 || relative_position(1,2) == 5)
437  {
438  return 0.0;
439  }
440 
441  for (int i = 0; i < 8; i++)
442  {
443  selected(i) = 0;
444  }
445 
446  // x = -1.0
447  r_pull_intersect(0,0) = -1.0;
448  r_pull_intersect(0,2) = vec1_pull(2) + \
449  (vec2_pull(2) - vec1_pull(2)) / (vec2_pull(0) - vec1_pull(0)) * (-1.0 - vec1_pull(0));
450  // x = 1.0
451  r_pull_intersect(1,0) = 1.0;
452  r_pull_intersect(1,2) = vec1_pull(2) + \
453  (vec2_pull(2) - vec1_pull(2)) / (vec2_pull(0) - vec1_pull(0)) * (1.0 - vec1_pull(0));
454  // z = -1.0
455  r_pull_intersect(4,0) = vec1_pull(0) + \
456  (vec2_pull(0) - vec1_pull(0)) / (vec2_pull(2) - vec1_pull(2)) * (-1.0 - vec1_pull(2));
457  r_pull_intersect(4,2) = -1.0;
458  // z = 1.0
459  r_pull_intersect(5,0) = vec1_pull(0) + \
460  (vec2_pull(0) - vec1_pull(0)) / (vec2_pull(2) - vec1_pull(2)) * (1.0 - vec1_pull(2));
461  r_pull_intersect(5,2) = 1.0;
462 
463  // vec1 itself
464  r_pull_intersect(6,0) = vec1_pull(0);
465  r_pull_intersect(6,1) = vec1_pull(1);
466  r_pull_intersect(6,2) = vec1_pull(2);
467 
468  // vec2 itself
469  r_pull_intersect(7,0) = vec2_pull(0);
470  r_pull_intersect(7,1) = vec2_pull(1);
471  r_pull_intersect(7,2) = vec2_pull(2);
472 
473  // if the point is inside, then it is the point selected.
474  if ((relative_position(0,0) == 2 || relative_position(0,0) == 3 || relative_position(0,0) == 4) && \
475  (relative_position(2,0) == 2 || relative_position(2,0) == 3 || relative_position(2,0) == 4))
476  {
477  selected(6) = 1;
478  }
479 
480  if ((relative_position(0,1) == 2 || relative_position(0,1) == 3 || relative_position(0,1) == 4) && \
481  (relative_position(2,1) == 2 || relative_position(2,1) == 3 || relative_position(2,1) == 4))
482  {
483  selected(7) = 1;
484  }
485 
486  // must be inside cubes to be selected
487  for (int i = 0; i <= 1; i++)
488  {
489  if (((r_pull_intersect(i,0) > vec1_pull(0) && r_pull_intersect(i,0) < vec2_pull(0)) || \
490  (r_pull_intersect(i,0) > vec2_pull(0) && r_pull_intersect(i,0) < vec1_pull(0))) && \
491  ((r_pull_intersect(i,2) > vec1_pull(2) && r_pull_intersect(i,2) < vec2_pull(2)) || \
492  (r_pull_intersect(i,2) > vec2_pull(2) && r_pull_intersect(i,2) < vec1_pull(2))) && \
493  (r_pull_intersect(i,2) >= -1.0 - epsilon && r_pull_intersect(i,2) <= 1.0 + epsilon))
494  {
495  selected(i) = 1;
496  }
497  }
498  for (int i = 4; i <= 5; i++)
499  {
500  if (((r_pull_intersect(i,0) > vec1_pull(0) && r_pull_intersect(i,0) < vec2_pull(0)) || \
501  (r_pull_intersect(i,0) > vec2_pull(0) && r_pull_intersect(i,0) < vec1_pull(0))) && \
502  ((r_pull_intersect(i,2) > vec1_pull(2) && r_pull_intersect(i,2) < vec2_pull(2)) || \
503  (r_pull_intersect(i,2) > vec2_pull(2) && r_pull_intersect(i,2) < vec1_pull(2))) && \
504  (r_pull_intersect(i,0) >= -1.0 - epsilon && r_pull_intersect(i,0) <= 1.0 + epsilon))
505  {
506  selected(i) = 1;
507  }
508  }
509 
510  selected_count = 0;
511  for (int i = 0; i < 8; i++)
512  {
513  if (selected(i) && selected_count == 0)
514  {
515  vec1_pull_seg(0) = r_pull_intersect(i,0);
516  vec1_pull_seg(2) = r_pull_intersect(i,2);
517  selected_count = 1;
518  continue;
519  }
520  if (selected(i) && selected_count == 1)
521  {
522  if (fabs(vec1_pull_seg(0) - r_pull_intersect(i,0)) < epsilon && \
523  fabs(vec1_pull_seg(2) - r_pull_intersect(i,2)) < epsilon)
524  {
525  continue;
526  }
527  else
528  {
529  vec2_pull_seg(0) = r_pull_intersect(i,0);
530  vec2_pull_seg(2) = r_pull_intersect(i,2);
531  selected_count = 2;
532  continue;
533  }
534  }
535  if (selected(i) && selected_count == 2)
536  {
537  if ((fabs(vec1_pull_seg(0) - r_pull_intersect(i,0)) < epsilon && \
538  fabs(vec1_pull_seg(2) - r_pull_intersect(i,2)) < epsilon) || \
539  (fabs(vec2_pull_seg(0) - r_pull_intersect(i,0)) < epsilon && \
540  fabs(vec2_pull_seg(2) - r_pull_intersect(i,2)) < epsilon))
541  {
542  continue;
543  }
544  else
545  {
546  MY_ERROR("In LDAD bond_function xz, More than 2 different points \
547  with different positions are being selected. \
548  This means a line intersecting with a cube has 3 points. \
549  Something wrong with the algorithm. Need to print and debug!");
550  }
551  }
552  }
553 
554  // No intersection .or. one point on the parallelepiped but another point is outside
555  if (selected_count == 0 || selected_count == 1)
556  {
557  return 0.0;
558  }
559  }
560  // degenerate to 2D yz
561  else if (degenerate(0) && !degenerate(1) && !degenerate(2))
562  {
563  relative_position(0,2) = PointLineRelationship((vec1_pull(0) + vec2_pull(0)) / 2.0);
564  if (relative_position(0,2) == 1 || relative_position(0,2) == 5)
565  {
566 
567  return 0.0;
568  }
569 
570  for (int i = 0; i < 8; i++)
571  {
572  selected(i) = 0;
573  }
574 
575  // y = -1.0
576  r_pull_intersect(2,1) = -1.0;
577  r_pull_intersect(2,2) = vec1_pull(2) + \
578  (vec2_pull(2) - vec1_pull(2)) / (vec2_pull(1) - vec1_pull(1)) * (-1.0 - vec1_pull(1));
579  // y = 1.0
580  r_pull_intersect(3,1) = 1.0;
581  r_pull_intersect(3,2) = vec1_pull(2) + \
582  (vec2_pull(2) - vec1_pull(2)) / (vec2_pull(1) - vec1_pull(1)) * (1.0 - vec1_pull(1));
583  // z = -1.0
584  r_pull_intersect(4,1) = vec1_pull(1) + \
585  (vec2_pull(1) - vec1_pull(1)) / (vec2_pull(2) - vec1_pull(2)) * (-1.0 - vec1_pull(2));
586  r_pull_intersect(4,2) = -1.0;
587  // z = 1.0
588  r_pull_intersect(5,1) = vec1_pull(1) + \
589  (vec2_pull(1) - vec1_pull(1)) / (vec2_pull(2) - vec1_pull(2)) * (1.0 - vec1_pull(2));
590  r_pull_intersect(5,2) = 1.0;
591 
592  // vec1 itself
593  r_pull_intersect(6,0) = vec1_pull(0);
594  r_pull_intersect(6,1) = vec1_pull(1);
595  r_pull_intersect(6,2) = vec1_pull(2);
596 
597  // vec2 itself
598  r_pull_intersect(7,0) = vec2_pull(0);
599  r_pull_intersect(7,1) = vec2_pull(1);
600  r_pull_intersect(7,2) = vec2_pull(2);
601 
602  // if the point is inside, then it is the point selected.
603  if ((relative_position(1,0) == 2 || relative_position(1,0) == 3 || relative_position(1,0) == 4) && \
604  (relative_position(2,0) == 2 || relative_position(2,0) == 3 || relative_position(2,0) == 4))
605  {
606  selected(6) = 1;
607  }
608 
609  if ((relative_position(1,1) == 2 || relative_position(1,1) == 3 || relative_position(1,1) == 4) && \
610  (relative_position(2,1) == 2 || relative_position(2,1) == 3 || relative_position(2,1) == 4))
611  {
612  selected(7) = 1;
613  }
614 
615  // must be inside cubes to be selected
616  for (int i = 2; i <= 3; i++)
617  {
618  if (((r_pull_intersect(i,1) > vec1_pull(1) && r_pull_intersect(i,1) < vec2_pull(1)) || \
619  (r_pull_intersect(i,1) > vec2_pull(1) && r_pull_intersect(i,1) < vec1_pull(1))) && \
620  ((r_pull_intersect(i,2) > vec1_pull(2) && r_pull_intersect(i,2) < vec2_pull(2)) || \
621  (r_pull_intersect(i,2) > vec2_pull(2) && r_pull_intersect(i,2) < vec1_pull(2))) && \
622  (r_pull_intersect(i,2) >= -1.0 - epsilon && r_pull_intersect(i,2) <= 1.0 + epsilon))
623  {
624  selected(i) = 1;
625  }
626  }
627  for (int i = 4; i <= 5; i++)
628  {
629  if (((r_pull_intersect(i,1) > vec1_pull(1) && r_pull_intersect(i,1) < vec2_pull(1)) || \
630  (r_pull_intersect(i,1) > vec2_pull(1) && r_pull_intersect(i,1) < vec1_pull(1))) && \
631  ((r_pull_intersect(i,2) > vec1_pull(2) && r_pull_intersect(i,2) < vec2_pull(2)) || \
632  (r_pull_intersect(i,2) > vec2_pull(2) && r_pull_intersect(i,2) < vec1_pull(2))) && \
633  (r_pull_intersect(i,1) >= -1.0 - epsilon && r_pull_intersect(i,1) <= 1.0 + epsilon))
634  {
635  selected(i) = 1;
636  }
637  }
638 
639  selected_count = 0;
640  for (int i = 0; i < 8; i++)
641  {
642  if (selected(i) && selected_count == 0)
643  {
644  vec1_pull_seg(1) = r_pull_intersect(i,1);
645  vec1_pull_seg(2) = r_pull_intersect(i,2);
646  selected_count = 1;
647  continue;
648  }
649  if (selected(i) && selected_count == 1)
650  {
651  if (fabs(vec1_pull_seg(1) - r_pull_intersect(i,1)) < epsilon && \
652  fabs(vec1_pull_seg(2) - r_pull_intersect(i,2)) < epsilon)
653  {
654  continue;
655  }
656  else
657  {
658  vec2_pull_seg(1) = r_pull_intersect(i,1);
659  vec2_pull_seg(2) = r_pull_intersect(i,2);
660  selected_count = 2;
661  continue;
662  }
663  }
664  if (selected(i) && selected_count == 2)
665  {
666  if ((fabs(vec1_pull_seg(1) - r_pull_intersect(i,1)) < epsilon && \
667  fabs(vec1_pull_seg(2) - r_pull_intersect(i,2)) < epsilon) || \
668  (fabs(vec2_pull_seg(1) - r_pull_intersect(i,1)) < epsilon && \
669  fabs(vec2_pull_seg(2) - r_pull_intersect(i,2)) < epsilon))
670  {
671  continue;
672  }
673  else
674  {
675  MY_ERROR("In LDAD bond_function yz, More than 2 different points \
676  with different positions are being selected. \
677  This means a line intersecting with a cube has 3 points. \
678  Something wrong with the algorithm. Need to print and debug!");
679  }
680  }
681  }
682 
683  // No intersection .or. one point on the parallelepiped but another point is outside
684  if (selected_count == 0 || selected_count == 1)
685  {
686 
687  return 0.0;
688  }
689  }
690  // 3D xyz
691  else if (!degenerate(0) && !degenerate(1) && !degenerate(2))
692  {
693  for (int i = 0; i < 8; i++)
694  {
695  selected(i) = 0;
696  }
697 
698  // The 8 points
699  // x = -1.0
700  r_pull_intersect(0,0) = -1.0;
701  r_pull_intersect(0,1) = vec1_pull(1) + \
702  (vec2_pull(1) - vec1_pull(1)) / (vec2_pull(0) - vec1_pull(0)) * (-1.0 - vec1_pull(0));
703  r_pull_intersect(0,2) = vec1_pull(2) + \
704  (vec2_pull(2) - vec1_pull(2)) / (vec2_pull(0) - vec1_pull(0)) * (-1.0 - vec1_pull(0));
705  // x = 1.0
706  r_pull_intersect(1,0) = 1.0;
707  r_pull_intersect(1,1) = vec1_pull(1) + \
708  (vec2_pull(1) - vec1_pull(1)) / (vec2_pull(0) - vec1_pull(0)) * (1.0 - vec1_pull(0));
709  r_pull_intersect(1,2) = vec1_pull(2) + \
710  (vec2_pull(2) - vec1_pull(2)) / (vec2_pull(0) - vec1_pull(0)) * (1.0 - vec1_pull(0));
711 
712  // y = -1.0
713  r_pull_intersect(2,0) = vec1_pull(0) + \
714  (vec2_pull(0) - vec1_pull(0)) / (vec2_pull(1) - vec1_pull(1)) * (-1.0 - vec1_pull(1));
715  r_pull_intersect(2,1) = -1.0;
716  r_pull_intersect(2,2) = vec1_pull(2) + \
717  (vec2_pull(2) - vec1_pull(2)) / (vec2_pull(1) - vec1_pull(1)) * (-1.0 - vec1_pull(1));
718 
719  // y = 1.0
720  r_pull_intersect(3,0) = vec1_pull(0) + \
721  (vec2_pull(0) - vec1_pull(0)) / (vec2_pull(1) - vec1_pull(1)) * (1.0 - vec1_pull(1));
722  r_pull_intersect(3,1) = 1.0;
723  r_pull_intersect(3,2) = vec1_pull(2) + \
724  (vec2_pull(2) - vec1_pull(2)) / (vec2_pull(1) - vec1_pull(1)) * (1.0 - vec1_pull(1));
725 
726  // z = -1.0
727  r_pull_intersect(4,0) = vec1_pull(0) + \
728  (vec2_pull(0) - vec1_pull(0)) / (vec2_pull(2) - vec1_pull(2)) * (-1.0 - vec1_pull(2));
729  r_pull_intersect(4,1) = vec1_pull(1) + \
730  (vec2_pull(1) - vec1_pull(1)) / (vec2_pull(2) - vec1_pull(2)) * (-1.0 - vec1_pull(2));
731  r_pull_intersect(4,2) = -1.0;
732 
733  // z = 1.0
734  r_pull_intersect(5,0) = vec1_pull(0) + \
735  (vec2_pull(0) - vec1_pull(0)) / (vec2_pull(2) - vec1_pull(2)) * (1.0 - vec1_pull(2));
736  r_pull_intersect(5,1) = vec1_pull(1) + \
737  (vec2_pull(1) - vec1_pull(1)) / (vec2_pull(2) - vec1_pull(2)) * (1.0 - vec1_pull(2));
738  r_pull_intersect(5,2) = 1.0;
739 
740  // vec1 itself
741  r_pull_intersect(6,0) = vec1_pull(0);
742  r_pull_intersect(6,1) = vec1_pull(1);
743  r_pull_intersect(6,2) = vec1_pull(2);
744 
745  // vec2 itself
746  r_pull_intersect(7,0) = vec2_pull(0);
747  r_pull_intersect(7,1) = vec2_pull(1);
748  r_pull_intersect(7,2) = vec2_pull(2);
749 
750  // if the point is inside, then it is the point selected.
751  if ((relative_position(0,0) == 2 || relative_position(0,0) == 3 || relative_position(0,0) == 4) && \
752  (relative_position(1,0) == 2 || relative_position(1,0) == 3 || relative_position(1,0) == 4) && \
753  (relative_position(2,0) == 2 || relative_position(2,0) == 3 || relative_position(2,0) == 4))
754  {
755  selected(6) = 1;
756  }
757 
758  if ((relative_position(0,1) == 2 || relative_position(0,1) == 3 || relative_position(0,1) == 4) && \
759  (relative_position(1,1) == 2 || relative_position(1,1) == 3 || relative_position(1,1) == 4) && \
760  (relative_position(2,1) == 2 || relative_position(2,1) == 3 || relative_position(2,1) == 4))
761  {
762  selected(7) = 1;
763  }
764 
765  for (int i = 0; i <= 5; i++)
766  {
767  // see whether the intersection point is between vec1 and vec2
768  if (((r_pull_intersect(i,0) > vec1_pull(0) && r_pull_intersect(i,0) < vec2_pull(0)) || \
769  (r_pull_intersect(i,0) > vec2_pull(0) && r_pull_intersect(i,0) < vec1_pull(0))) && \
770  ((r_pull_intersect(i,1) > vec1_pull(1) && r_pull_intersect(i,1) < vec2_pull(1)) || \
771  (r_pull_intersect(i,1) > vec2_pull(1) && r_pull_intersect(i,1) < vec1_pull(1))) && \
772  ((r_pull_intersect(i,2) > vec1_pull(2) && r_pull_intersect(i,2) < vec2_pull(2)) || \
773  (r_pull_intersect(i,2) > vec2_pull(2) && r_pull_intersect(i,2) < vec1_pull(2))) && \
774  (r_pull_intersect(i,0) >= -1.0 - epsilon && r_pull_intersect(i,0) <= 1.0 + epsilon) && \
775  (r_pull_intersect(i,1) >= -1.0 - epsilon && r_pull_intersect(i,1) <= 1.0 + epsilon) && \
776  (r_pull_intersect(i,2) >= -1.0 - epsilon && r_pull_intersect(i,2) <= 1.0 + epsilon))
777  {
778  selected(i) = 1;
779  }
780  }
781 
782  selected_count = 0;
783 
784  for (int i = 0; i < 8; i++)
785  {
786  if (selected(i) && selected_count == 0)
787  {
788  vec1_pull_seg(0) = r_pull_intersect(i,0);
789  vec1_pull_seg(1) = r_pull_intersect(i,1);
790  vec1_pull_seg(2) = r_pull_intersect(i,2);
791  selected_count = 1;
792  continue;
793  }
794 
795  if (selected(i) && selected_count == 1)
796  {
797  if (fabs(vec1_pull_seg(0) - r_pull_intersect(i,0)) < epsilon && \
798  fabs(vec1_pull_seg(1) - r_pull_intersect(i,1)) < epsilon && \
799  fabs(vec1_pull_seg(2) - r_pull_intersect(i,2)) < epsilon)
800  {
801  continue;
802  }
803  else
804  {
805  vec2_pull_seg(0) = r_pull_intersect(i,0);
806  vec2_pull_seg(1) = r_pull_intersect(i,1);
807  vec2_pull_seg(2) = r_pull_intersect(i,2);
808  selected_count = 2;
809  continue;
810  }
811  }
812 
813  if (selected(i) && selected_count == 2)
814  {
815  if ((fabs(vec1_pull_seg(0) - r_pull_intersect(i,0)) < epsilon && \
816  fabs(vec1_pull_seg(1) - r_pull_intersect(i,1)) < epsilon && \
817  fabs(vec1_pull_seg(2) - r_pull_intersect(i,2)) < epsilon) || \
818  (fabs(vec2_pull_seg(0) - r_pull_intersect(i,0)) < epsilon && \
819  fabs(vec2_pull_seg(1) - r_pull_intersect(i,1)) < epsilon && \
820  fabs(vec2_pull_seg(2) - r_pull_intersect(i,2)) < epsilon))
821  {
822  continue;
823  }
824  else
825  {
826  MY_ERROR("In LDAD bond_function xyz, More than 2 different points \
827  with different positions are being selected. \
828  This means a line intersecting with a cube has 3 points. \
829  Something wrong with the algorithm. Need to print and debug!");
830  }
831  }
832  }
833 
834  // No intersection .or. one point on the parallelepiped but another point is outside
835  if (selected_count == 0 || selected_count == 1)
836  {
837  return 0.0;
838  }
839  }
840 
841  vec1_push_seg = ldadVectors * vec1_pull_seg.transpose();
842  vec2_push_seg = ldadVectors * vec2_pull_seg.transpose();
843  vec12_push_seg = vec2_push_seg - vec1_push_seg;
844  total_length = vec12_push_seg.norm();
845 
846  // use oneDFunction.integrate(vec1_pull_seg, vec2_pull_seg) -> helper function;
847  return total_length * oneDFunction.integrate(vec1_pull_seg, vec2_pull_seg) * normalizer / distance;
848 }
849 
850 int PointLineRelationship(const double& p)
851 {
852  if (p < -1.0 - epsilon)
853  {
854  return 1;
855  }
856  else if (p >= -1.0 - epsilon && p <= -1.0 + epsilon)
857  {
858  return 2;
859  }
860  else if (p > -1.0 + epsilon && p < 1.0 - epsilon)
861  {
862  return 3;
863  }
864  else if (p >= 1.0 - epsilon && p <= 1.0 + epsilon)
865  {
866  return 4;
867  }
868  else if (p > 1.0 + epsilon)
869  {
870  return 5;
871  }
872  else
873  {
874  std::cout << "Point " << p << std::endl;
875  MY_ERROR("In PointLineRelationship, the above point did not fall into any range.");
876  }
877 }
epsilon
const double epsilon
Definition: typedef.h:73
Method< MethodLdad< T > >::averagingDomainSize
double averagingDomainSize
Definition: Method.h:24
MethodLdad::operator()
double operator()(const Vector3d &vec) const
Definition: MethodLdad.cpp:46
MethodLdad::~MethodLdad
virtual ~MethodLdad()
Definition: MethodLdad.cpp:41
MethodLdad::bondFunction
double bondFunction(const Vector3d &vec1, const Vector3d &vec2) const
Definition: MethodLdad.cpp:55
PointLineRelationship
int PointLineRelationship(const double &p)
Definition: MethodLdad.cpp:850
MY_ERROR
#define MY_ERROR(message)
Definition: typedef.h:17
norm
double norm(double const *a)
Definition: helper.hpp:20
Vector3d
Eigen::Matrix< double, 1, DIM, Eigen::RowMajor > Vector3d
Definition: typedef.h:60
Matrix3d
Eigen::Matrix< double, DIM, DIM, Eigen::RowMajor > Matrix3d
Definition: typedef.h:56
MethodLdad::MethodLdad
MethodLdad(const Matrix3d &ldadVectors)
Definition: MethodLdad.cpp:16
Matrix3i
Eigen::Matrix< int, DIM, DIM, Eigen::RowMajor > Matrix3i
Definition: typedef.h:57
Vector3i
Eigen::Matrix< int, 1, DIM, Eigen::RowMajor > Vector3i
Definition: typedef.h:61
ArrayXXd
Eigen::Array< double, Eigen::Dynamic, Eigen::Dynamic > ArrayXXd
Definition: typedef.h:62
VectorXi
Eigen::Matrix< int, 1, Eigen::Dynamic, Eigen::RowMajor > VectorXi
Definition: typedef.h:58