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