Chaste Commit::baa90ac2819b962188b7562f2326be23c47859a7
ObsoleteBoxCollection.cpp
1/*
2
3Copyright (c) 2005-2024, University of Oxford.
4All rights reserved.
5
6University of Oxford means the Chancellor, Masters and Scholars of the
7University of Oxford, having an administrative office at Wellington
8Square, Oxford OX1 2JD, UK.
9
10This file is part of Chaste.
11
12Redistribution and use in source and binary forms, with or without
13modification, are permitted provided that the following conditions are met:
14 * Redistributions of source code must retain the above copyright notice,
15 this list of conditions and the following disclaimer.
16 * Redistributions in binary form must reproduce the above copyright notice,
17 this list of conditions and the following disclaimer in the documentation
18 and/or other materials provided with the distribution.
19 * Neither the name of the University of Oxford nor the names of its
20 contributors may be used to endorse or promote products derived from this
21 software without specific prior written permission.
22
23THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
24AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
25IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
26ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
27LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
28CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE
29GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
30HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
31LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT
32OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
33
34*/
35
36#include "ObsoleteBoxCollection.hpp"
37#include "Exception.hpp"
38
40// ObsoleteBoxCollection methods
42
43// Static member for "fudge factor" is instantiated here
44template<unsigned DIM>
45const double ObsoleteBoxCollection<DIM>::msFudge = 5e-14;
46
47template<unsigned DIM>
48ObsoleteBoxCollection<DIM>::ObsoleteBoxCollection(double boxWidth, c_vector<double, 2 * DIM> domainSize, bool isPeriodicInX,
49 bool isPeriodicInY, bool isPeriodicInZ)
50 : mDomainSize(domainSize),
51 mBoxWidth(boxWidth)
52{
53 // Populate mIsDomainPeriodic
54 switch (DIM)
55 {
56 case 1:
57 {
58 mIsDomainPeriodic[0] = isPeriodicInX;
59 break;
60 }
61 case 2:
62 {
63 mIsDomainPeriodic[0] = isPeriodicInX;
64 mIsDomainPeriodic[1] = isPeriodicInY;
65 break;
66 }
67 case 3:
68 {
69 mIsDomainPeriodic[0] = isPeriodicInX;
70 mIsDomainPeriodic[1] = isPeriodicInY;
71 mIsDomainPeriodic[2] = isPeriodicInZ;
72 break;
73 }
74 default:
75 {
77 }
78 }
79
80 // Calculating the number of boxes in each direction and total number of boxes.
81 unsigned num_boxes = 1;
83 for (unsigned i = 0; i < DIM; i++)
84 {
86 mNumBoxesEachDirection(i) = (unsigned) floor((domainSize(2 * i + 1) - domainSize(2 * i)) / boxWidth + msFudge) + 1;
87 num_boxes *= mNumBoxesEachDirection(i);
88 }
89
90 // Create the correct number of boxes
91 mBoxes.resize(num_boxes);
92}
93
94template<unsigned DIM>
96{
97 for (unsigned i = 0; i < mBoxes.size(); i++)
98 {
99 mBoxes[i].ClearNodes();
100 }
101}
102
103template<unsigned DIM>
105{
106 // Get the location of the node
107 c_vector<double, DIM> location = pNode->rGetLocation();
108 return CalculateContainingBox(location);
109}
110
111template<unsigned DIM>
112unsigned ObsoleteBoxCollection<DIM>::CalculateContainingBox(c_vector<double, DIM>& rLocation)
113{
114 // Confirm that the location is in the domain
115 for (unsigned i = 0; i < DIM; i++)
116 {
117 if (!(rLocation[i] >= mDomainSize[2*i] && rLocation[i] <= mDomainSize[2*i + 1]))
118 {
119 std::stringstream location_error;
120 location_error << "Location in dimension " << i << " is " << rLocation[i] << " which is not in the domain ["
121 << mDomainSize[2*i] << ", " << mDomainSize[2*i + 1] << "].";
122
123 EXCEPTION(location_error.str());
124 }
125 }
126
127 // Compute the containing box index in each dimension
128 c_vector<int, DIM> containing_box_indices;
129 for (unsigned i = 0; i < DIM; i++)
130 {
131 containing_box_indices[i] = (int) floor((rLocation[i] - mDomainSize(2 * i) + msFudge) / mBoxWidth);
132 }
133
134 return GetLinearIndex(containing_box_indices);
135}
136
137template<unsigned DIM>
139{
140 assert(boxIndex < mBoxes.size());
141 return mBoxes[boxIndex];
142}
144template<unsigned DIM>
146{
147 return mBoxes.size();
148}
150template<unsigned DIM>
151unsigned ObsoleteBoxCollection<DIM>::GetLinearIndex(c_vector<int, DIM> gridIndices)
152{
153 /*
154 * This function may be passed values outside the range in one or more
155 * dimensions in the case of a periodic domain. We therefore assume that
156 * these values represent a situation with periodicity and adjust them
157 * accordingly before calculating the linear index.
158 */
159
160 // Adjust for periodicity if necessary
161 for (unsigned dim = 0; dim < DIM; dim++)
162 {
163 // Check for values too large
164 if (gridIndices(dim) >= (int)mNumBoxesEachDirection(dim))
165 {
166 assert(mIsDomainPeriodic(dim));
167 gridIndices(dim) -= (int)mNumBoxesEachDirection(dim);
168 }
169
170 // Check for negative values
171 else if (gridIndices(dim) < 0)
172 {
173 assert(mIsDomainPeriodic(dim));
174 gridIndices(dim) += (int)mNumBoxesEachDirection(dim);
175 }
176 }
177
178 // Calculate linear index
179 unsigned linear_index;
180
181 switch (DIM)
182 {
183 case 1:
184 {
185 linear_index = gridIndices(0);
186 break;
187 }
188 case 2:
189 {
190 linear_index = gridIndices(0) +
191 gridIndices(1) * mNumBoxesEachDirection(0);
192 break;
193 }
194 case 3:
195 {
196 linear_index = gridIndices(0) +
197 gridIndices(1) * mNumBoxesEachDirection(0) +
198 gridIndices(2) * mNumBoxesEachDirection(0) * mNumBoxesEachDirection(1);
199 break;
200 }
201 default:
202 {
204 }
205 }
206
207 return linear_index;
208}
209
210template<unsigned DIM>
211c_vector<unsigned, DIM> ObsoleteBoxCollection<DIM>::GetGridIndices(unsigned linearIndex)
212{
213 c_vector<unsigned, DIM> grid_indices {};
214
215 if constexpr (DIM == 1)
216 {
217 grid_indices(0) = linearIndex;
218 }
219 else if constexpr (DIM == 2)
220 {
221 const unsigned num_x = mNumBoxesEachDirection(0);
222 grid_indices(0) = linearIndex % num_x;
223 grid_indices(1) = (linearIndex - grid_indices(0)) / num_x;
224 }
225 else if constexpr (DIM == 3)
226 {
227 const unsigned num_x = mNumBoxesEachDirection(0);
228 const unsigned num_xy = mNumBoxesEachDirection(0)*mNumBoxesEachDirection(1);
229 grid_indices(0) = linearIndex % num_x;
230 grid_indices(1) = (linearIndex % num_xy - grid_indices(0)) / num_x;
231 grid_indices(2) = linearIndex / num_xy;
232 }
233
234 return grid_indices;
235}
236
237template<unsigned DIM>
238bool ObsoleteBoxCollection<DIM>::IsBoxInDomain(c_vector<unsigned, DIM> gridIndices)
239{
240 /*
241 * We assume that, for a given dimension, any location is in the domain
242 * if the domain is periodic in that dimension.
243 *
244 * This method can only return false in the case of one or more
245 * dimensions being aperiodic and a location in one of those dimensions
246 * being outside the range.
247 */
248
249 bool is_in_domain = true;
250
251 for (unsigned dim = 0; dim < DIM; dim++)
252 {
253 if (!mIsDomainPeriodic(dim))
254 {
255 if (gridIndices(dim) >= mNumBoxesEachDirection(dim))
256 {
257 is_in_domain = false;
258 }
259 }
260 }
261
262 return is_in_domain;
263}
264
265template<unsigned DIM>
266c_vector<bool,DIM> ObsoleteBoxCollection<DIM>::IsIndexPenultimate(c_vector<unsigned, DIM> gridIndices)
267{
268 c_vector<bool,DIM> is_penultimate;
269
270 for (unsigned dim = 0 ; dim < DIM ; dim++)
271 {
272 is_penultimate(dim) = (gridIndices(dim) == mNumBoxesEachDirection(dim) - 2);
273 }
274
275 return is_penultimate;
276}
277
278template<unsigned DIM>
280{
281 // Populate a list of half the neighbours in this number of dimensions
282 std::vector<c_vector<int,DIM> > neighbours;
283
284 switch (DIM)
285 {
286 case 1:
287 {
288 // Just one neighbour, plus the current box (zero vector)
289 neighbours = std::vector<c_vector<int,DIM> >(2);
290
291 neighbours[0](0) = 0; // current box
292 neighbours[1](0) = 1; // right
293
294 break;
295 }
296 case 2:
297 {
298 // Four neighbours, plus the current box (zero vector)
299 neighbours = std::vector<c_vector<int,DIM> >(5);
300
301 neighbours[0](0) = 0; neighbours[0](1) = 1; // up
302 neighbours[1](0) = 1; neighbours[1](1) = 1; // up right
303 neighbours[2](0) = 1; neighbours[2](1) = 0; // right
304 neighbours[3](0) = 1; neighbours[3](1) =-1; // down right
305 neighbours[4](0) = 0; neighbours[4](1) = 0; // current box
306
307 break;
308 }
309 case 3:
310 {
311 /*
312 * We need to pick 13 neighbours in such a way as to ensure all interactions are captured,
313 * in addition to the current box.
314 *
315 * The 26 cubes on the outside of a 3x3 arrangement either adjacent to a whole face of
316 * the central cube, only an edge of the central cube, or only a vertex of the central
317 * cube:
318 *
319 * 6 contact faces, and come in the following 3 pairs:
320 * +x -x (1, 0, 0)
321 * +y -y (0, 1, 0)
322 * +z -z (0, 0, 1)
323 *
324 * 12 contact edges, and come in the following 6 pairs:
325 * +x+y -x-y (1, 1, 0)
326 * +x+z -x-z (1, 0, 1)
327 * +x-y -x+y (1,-1, 0)
328 * +x-z -x+z (1, 0,-1)
329 * +y+z -y-z (0, 1, 1)
330 * +y-z -y+z (0, 1,-1)
331 *
332 * 8 contact vertices, and come in the following 4 pairs:
333 * +x+y+z -x-y-z (1, 1, 1)
334 * +x+y-z -x-y+z (1, 1,-1)
335 * +x-y+z -x+y-z (1,-1, 1)
336 * +x-y-z -x+y+z (1,-1,-1)
337 *
338 * We will simply pick the 13 from the left-hand column above, which can be represented
339 * as an integer offset in three dimensions from the middle cube, as shown in the third
340 * column above.
341 */
342
343 neighbours = std::vector<c_vector<int,DIM> >(14);
344
345 neighbours[0](0) = 1; neighbours[0](1) = 0; neighbours[0](2) = 0;
346 neighbours[1](0) = 0; neighbours[1](1) = 1; neighbours[1](2) = 0;
347 neighbours[2](0) = 0; neighbours[2](1) = 0; neighbours[2](2) = 1;
348 neighbours[3](0) = 1; neighbours[3](1) = 1; neighbours[3](2) = 0;
349 neighbours[4](0) = 1; neighbours[4](1) = 0; neighbours[4](2) = 1;
350 neighbours[5](0) = 1; neighbours[5](1) =-1; neighbours[5](2) = 0;
351 neighbours[6](0) = 1; neighbours[6](1) = 0; neighbours[6](2) =-1;
352 neighbours[7](0) = 0; neighbours[7](1) = 1; neighbours[7](2) = 1;
353 neighbours[8](0) = 0; neighbours[8](1) = 1; neighbours[8](2) =-1;
354 neighbours[9](0) = 1; neighbours[9](1) = 1; neighbours[9](2) = 1;
355 neighbours[10](0) = 1; neighbours[10](1) = 1; neighbours[10](2) =-1;
356 neighbours[11](0) = 1; neighbours[11](1) =-1; neighbours[11](2) = 1;
357 neighbours[12](0) = 1; neighbours[12](1) =-1; neighbours[12](2) =-1;
358
359 neighbours[13](0) = 0; neighbours[13](1) = 0; neighbours[13](2) = 0; // current box
360
361 break;
362 }
363 default:
364 {
366 }
367 }
368
369 // Pass the list of possible neighbours to SetupLocalBoxes()
370 SetupLocalBoxes(neighbours);
371}
372
373template<unsigned DIM>
375{
376 // Populate a list of all neighbours in this number of dimensions
377 std::vector<c_vector<int, DIM> > neighbours;
378
379 switch (DIM)
380 {
381 case 1:
382 {
383 // 3 neighbours
384 neighbours.clear();
385
386 for (int x_dim = -1 ; x_dim < 2 ; x_dim++)
387 {
388 c_vector<int, DIM> neighbour;
389 neighbour(0) = x_dim;
390
391 neighbours.push_back(neighbour);
392 }
393
394 break;
395 }
396 case 2:
397 {
398 // 3x3 neighbours
399 neighbours.clear();
400
401 for (int x_dim = -1 ; x_dim < 2 ; x_dim++)
402 {
403 for (int y_dim = -1 ; y_dim < 2 ; y_dim++)
404 {
405 c_vector<int, DIM> neighbour;
406 neighbour(0) = x_dim;
407 neighbour(1) = y_dim;
408
409 neighbours.push_back(neighbour);
410 }
411 }
412
413 break;
414 }
415 case 3:
416 {
417 // 3x3x3 neighbours
418 neighbours.clear();
419
420 for (int x_dim = -1 ; x_dim < 2 ; x_dim++)
421 {
422 for (int y_dim = -1 ; y_dim < 2 ; y_dim++)
423 {
424 for (int z_dim = -1 ; z_dim < 2 ; z_dim++)
425 {
426 c_vector<int, DIM> neighbour;
427 neighbour(0) = x_dim;
428 neighbour(1) = y_dim;
429 neighbour(2) = z_dim;
430
431 neighbours.push_back(neighbour);
432 }
433 }
434 }
435
436 break;
437 }
438 default:
439 {
441 }
442 }
443
444 // Pass the list of possible neighbours to SetupLocalBoxes()
445 SetupLocalBoxes(neighbours);
446}
447
448template<unsigned DIM>
449void ObsoleteBoxCollection<DIM>::SetupLocalBoxes(const std::vector<c_vector<int, DIM> >& rNeighbours)
450{
451 mLocalBoxes.clear();
452
453 // Loop over all boxes and add all necessary local boxes
454 for (unsigned box_idx = 0 ; box_idx < mBoxes.size() ; box_idx++)
455 {
456 std::set<unsigned> local_boxes;
457
458 // The grid indices (i), (i,j) or (i,j,k) depending on DIM, corresponding to box_idx
459 c_vector<unsigned, DIM> grid_indices = GetGridIndices(box_idx);
460
461 // Check all neighbours (1 or 2 in 1D, 4 or 8 in 2D, 13 or 26 in 3D) and add them as local
462 // boxes if they are in the domain
463 for (unsigned neighbour = 0 ; neighbour < rNeighbours.size() ; neighbour++)
464 {
465 if (IsBoxInDomain(grid_indices + rNeighbours[neighbour]))
466 {
467 local_boxes.insert(GetLinearIndex(grid_indices + rNeighbours[neighbour]));
468 }
469 }
470
471 /*
472 * If we are in a periodic domain, we need to add additional boxes if the
473 * current box is in the penultimate position of any dimension.
474 *
475 * The additional boxes we need to add are just all the neighbours of the box
476 * in the final position in that dimension. Because we use set::insert, it
477 * doesn't matter if we try and add the same box multiple times (only one
478 * copy of the index will be included).
479 */
480
481 // Find if the current indices are penultimate in any dimension in which collection is periodic
482 c_vector<bool, DIM> is_penultimate_periodic = IsIndexPenultimate(grid_indices);
483 unsigned num_penultimate_periodic_dimensions = 0;
484
485 for (unsigned dim = 0 ; dim < DIM ; dim++)
486 {
487 is_penultimate_periodic(dim) = is_penultimate_periodic(dim) && mIsDomainPeriodic(dim);
488 if (is_penultimate_periodic(dim))
489 {
490 num_penultimate_periodic_dimensions++;
491 }
492 }
493
494 // First, deal with at least one penultimate dimension
495 if (num_penultimate_periodic_dimensions > 0)
496 {
497 // Loop over each dimension. If penultimate in dimension DIM move one box and add all neighbours
498 for (unsigned dim = 0 ; dim < DIM ; dim++)
499 {
500 if (is_penultimate_periodic(dim))
501 {
502 // If we're penultimate in dimension dim, move to the final location and add each neighbour as before.
503 // The call to IsBoxInDomain() will handle periodicity.
504 // Note that we define this vector before setting it as otherwise the profiling build will break (see #2367)
505 c_vector<unsigned, DIM> ultimate_indices; // Initialisation on next line for profile compiler
506 ultimate_indices = grid_indices;
507 ultimate_indices(dim) ++;
508
509 for (unsigned neighbour = 0 ; neighbour < rNeighbours.size() ; neighbour++)
510 {
511 if (IsBoxInDomain(ultimate_indices + rNeighbours[neighbour]))
512 {
513 local_boxes.insert(GetLinearIndex(ultimate_indices + rNeighbours[neighbour]));
514 }
515 }
516 }
517 }
518 }
519
520 // The final consideration is cases of multiple dimensions of periodicity.
521 if (num_penultimate_periodic_dimensions > 1)
522 {
523 switch (DIM)
524 {
525 case 2:
526 {
527 // Must have x and y periodicity - just have to add one extra box
528 assert(mIsDomainPeriodic(0) && mIsDomainPeriodic(1));
529
530 local_boxes.insert(0);
531
532 break;
533 }
534 case 3:
535 {
536 // Could have X and Y, X and Z, Y and Z, or X, Y and Z periodicity
537
538 // X and Y
539 if (is_penultimate_periodic(0) && is_penultimate_periodic(1))
540 {
541 c_vector<unsigned, DIM> ultimate_indices;
542 ultimate_indices = grid_indices;
543 ultimate_indices(0) ++;
544 ultimate_indices(1) ++;
545
546 for (unsigned neighbour = 0 ; neighbour < rNeighbours.size() ; neighbour++)
547 {
548 if (IsBoxInDomain(ultimate_indices + rNeighbours[neighbour]))
549 {
550 local_boxes.insert(GetLinearIndex(ultimate_indices + rNeighbours[neighbour]));
551 }
552 }
553 }
554
555 // X and Z
556 if (is_penultimate_periodic(0) && is_penultimate_periodic(2))
557 {
558 c_vector<unsigned, DIM> ultimate_indices;
559 ultimate_indices = grid_indices;
560 ultimate_indices(0) ++;
561 ultimate_indices(2) ++;
562
563 for (unsigned neighbour = 0 ; neighbour < rNeighbours.size() ; neighbour++)
564 {
565 if (IsBoxInDomain(ultimate_indices + rNeighbours[neighbour]))
566 {
567 local_boxes.insert(GetLinearIndex(ultimate_indices + rNeighbours[neighbour]));
568 }
569 }
570 }
571
572 // Y and Z
573 if (is_penultimate_periodic(1) && is_penultimate_periodic(2))
574 {
575 c_vector<unsigned, DIM> ultimate_indices;
576 ultimate_indices = grid_indices;
577 ultimate_indices(1) ++;
578 ultimate_indices(2) ++;
579
580 for (unsigned neighbour = 0 ; neighbour < rNeighbours.size() ; neighbour++)
581 {
582 if (IsBoxInDomain(ultimate_indices + rNeighbours[neighbour]))
583 {
584 local_boxes.insert(GetLinearIndex(ultimate_indices + rNeighbours[neighbour]));
585 }
586 }
587 }
588
589 // X Y and Z
590 if (num_penultimate_periodic_dimensions == 3)
591 {
592 assert(mIsDomainPeriodic(0) && mIsDomainPeriodic(1) && mIsDomainPeriodic(1));
593 local_boxes.insert(0);
594 }
595
596 break;
597 }
598 default:
599 {
601 }
602 }
603 }
604
605 // Add the local boxes to the member vector
606 mLocalBoxes.push_back(local_boxes);
607 }
608}
609
610template<unsigned DIM>
611std::set<unsigned>& ObsoleteBoxCollection<DIM>::rGetLocalBoxes(unsigned boxIndex)
612{
613 assert(boxIndex < mLocalBoxes.size());
614 return mLocalBoxes[boxIndex];
615}
616
617template<unsigned DIM>
618const c_vector<double, 2 * DIM>& ObsoleteBoxCollection<DIM>::rGetDomainSize() const
619{
620 return mDomainSize;
621}
622
623template<unsigned DIM>
625 std::vector<std::pair<Node<DIM>*, Node<DIM>*> >& rNodePairs)
626{
627 rNodePairs.clear();
628
629 // Ensure all boxes are empty
630 EmptyBoxes();
631
632 // Create an empty set of neighbours for each node, and add each node to its correct box
633 for (unsigned node_index = 0; node_index < rNodes.size(); node_index++)
634 {
635 rNodes[node_index]->ClearNeighbours();
636
637 unsigned box_index = CalculateContainingBox(rNodes[node_index]);
638 mBoxes[box_index].AddNode(rNodes[node_index]);
639 }
640
641 for (unsigned i = 0; i < rNodes.size(); i++)
642 {
643 Node<DIM>* this_node = rNodes[i];
644 unsigned node_index = this_node->GetIndex();
645
646 // Get the box containing this node
647 unsigned this_node_box_index = CalculateContainingBox(this_node);
648
649 // Get the local boxes to this node
650 std::set<unsigned> local_boxes_indices = rGetLocalBoxes(this_node_box_index);
651
652 // Loop over all the local boxes
653 for (std::set<unsigned>::iterator box_iter = local_boxes_indices.begin(); box_iter != local_boxes_indices.end();
654 box_iter++)
655 {
656 // Get the set of nodes contained in this box
657 std::set<Node<DIM>*>& r_contained_nodes = mBoxes[*box_iter].rGetNodesContained();
658
659 // Loop over these nodes
660 for (typename std::set<Node<DIM>*>::iterator node_iter = r_contained_nodes.begin();
661 node_iter != r_contained_nodes.end(); ++node_iter)
662 {
663 // Get the index of the other node
664 unsigned other_node_index = (*node_iter)->GetIndex();
665
666 // If we're in the same box, then take care not to store the node pair twice
667 if (*box_iter == this_node_box_index)
668 {
669 if (other_node_index > this_node->GetIndex())
670 {
671 rNodePairs.push_back(std::pair<Node<DIM>*, Node<DIM>*>(this_node, (*node_iter)));
672 this_node->AddNeighbour(other_node_index);
673 (*node_iter)->AddNeighbour(node_index);
674 }
675 }
676 else
677 {
678 rNodePairs.push_back(std::pair<Node<DIM>*, Node<DIM>*>(this_node, (*node_iter)));
679 this_node->AddNeighbour(other_node_index);
680 (*node_iter)->AddNeighbour(node_index);
681 }
682 }
683 }
684 }
685
686 for (unsigned i = 0; i < rNodes.size(); i++)
687 {
688 rNodes[i]->RemoveDuplicateNeighbours();
689 }
690}
691
693
694template class ObsoleteBoxCollection<1>;
695template class ObsoleteBoxCollection<2>;
696template class ObsoleteBoxCollection<3>;
#define EXCEPTION(message)
#define NEVER_REACHED
Definition Box.hpp:51
Definition Node.hpp:59
void AddNeighbour(unsigned index)
Definition Node.cpp:328
const c_vector< double, SPACE_DIM > & rGetLocation() const
Definition Node.cpp:139
unsigned GetIndex() const
Definition Node.cpp:158
bool IsBoxInDomain(c_vector< unsigned, DIM > gridIndices)
std::vector< Box< DIM > > mBoxes
void SetupLocalBoxes(const std::vector< c_vector< int, DIM > > &rNeighbours)
c_vector< bool, DIM > mIsDomainPeriodic
Box< DIM > & rGetBox(unsigned boxIndex)
c_vector< unsigned, DIM > GetGridIndices(unsigned linearIndex)
void CalculateNodePairs(const std::vector< Node< DIM > * > &rNodes, std::vector< std::pair< Node< DIM > *, Node< DIM > * > > &rNodePairs)
c_vector< unsigned, DIM > mNumBoxesEachDirection
c_vector< bool, DIM > IsIndexPenultimate(c_vector< unsigned, DIM > gridIndices)
std::set< unsigned > & rGetLocalBoxes(unsigned boxIndex)
unsigned GetLinearIndex(c_vector< int, DIM > gridIndices)
ObsoleteBoxCollection(double boxWidth, c_vector< double, 2 *DIM > domainSize, bool isPeriodicInX=false, bool isPeriodicInY=false, bool isPeriodicInZ=false)
const c_vector< double, 2 *DIM > & rGetDomainSize() const
unsigned CalculateContainingBox(Node< DIM > *pNode)