MDStressLab++
Loading...
Searching...
No Matches
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
13int PointLineRelationship(const double& p);
14
15template<typename T>
16MethodLdad<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
40template<typename T>
41MethodLdad<T>::~MethodLdad() {
42 // TODO Auto-generated destructor stub
43}
44
45template<typename T>
46double 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
54template<typename T>
55double 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
850int 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}
PointLineRelationship
int PointLineRelationship(const double &p)
Definition MethodLdad.cpp:850
MethodLdad::MethodLdad
MethodLdad(const Matrix3d &ldadVectors)
Definition MethodLdad.cpp:16
MethodLdad::~MethodLdad
virtual ~MethodLdad()
Definition MethodLdad.cpp:41
MethodLdad::operator()
double operator()(const Vector3d &vec) const
Definition MethodLdad.cpp:46
MethodLdad::bondFunction
double bondFunction(const Vector3d &vec1, const Vector3d &vec2) const
Definition MethodLdad.cpp:55
norm
double norm(double const *a)
Definition helper.hpp:20
MY_ERROR
#define MY_ERROR(message)
Definition typedef.h:17
Matrix3d
Eigen::Matrix< double, DIM, DIM, Eigen::RowMajor > Matrix3d
Definition typedef.h:56
VectorXi
Eigen::Matrix< int, 1, Eigen::Dynamic, Eigen::RowMajor > VectorXi
Definition typedef.h:58
Vector3d
Eigen::Matrix< double, 1, DIM, Eigen::RowMajor > Vector3d
Definition typedef.h:60
Vector3i
Eigen::Matrix< int, 1, DIM, Eigen::RowMajor > Vector3i
Definition typedef.h:61
Matrix3i
Eigen::Matrix< int, DIM, DIM, Eigen::RowMajor > Matrix3i
Definition typedef.h:57
epsilon
const double epsilon
Definition typedef.h:73
ArrayXXd
Eigen::Array< double, Eigen::Dynamic, Eigen::Dynamic > ArrayXXd
Definition typedef.h:62