213 c_vector<unsigned, DIM> grid_indices {};
215 if constexpr (DIM == 1)
217 grid_indices(0) = linearIndex;
219 else if constexpr (DIM == 2)
221 const unsigned num_x = mNumBoxesEachDirection(0);
222 grid_indices(0) = linearIndex % num_x;
223 grid_indices(1) = (linearIndex - grid_indices(0)) / num_x;
225 else if constexpr (DIM == 3)
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;
282 std::vector<c_vector<int,DIM> > neighbours;
289 neighbours = std::vector<c_vector<int,DIM> >(2);
291 neighbours[0](0) = 0;
292 neighbours[1](0) = 1;
299 neighbours = std::vector<c_vector<int,DIM> >(5);
301 neighbours[0](0) = 0; neighbours[0](1) = 1;
302 neighbours[1](0) = 1; neighbours[1](1) = 1;
303 neighbours[2](0) = 1; neighbours[2](1) = 0;
304 neighbours[3](0) = 1; neighbours[3](1) =-1;
305 neighbours[4](0) = 0; neighbours[4](1) = 0;
343 neighbours = std::vector<c_vector<int,DIM> >(14);
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;
359 neighbours[13](0) = 0; neighbours[13](1) = 0; neighbours[13](2) = 0;
370 SetupLocalBoxes(neighbours);
377 std::vector<c_vector<int, DIM> > neighbours;
386 for (
int x_dim = -1 ; x_dim < 2 ; x_dim++)
388 c_vector<int, DIM> neighbour;
389 neighbour(0) = x_dim;
391 neighbours.push_back(neighbour);
401 for (
int x_dim = -1 ; x_dim < 2 ; x_dim++)
403 for (
int y_dim = -1 ; y_dim < 2 ; y_dim++)
405 c_vector<int, DIM> neighbour;
406 neighbour(0) = x_dim;
407 neighbour(1) = y_dim;
409 neighbours.push_back(neighbour);
420 for (
int x_dim = -1 ; x_dim < 2 ; x_dim++)
422 for (
int y_dim = -1 ; y_dim < 2 ; y_dim++)
424 for (
int z_dim = -1 ; z_dim < 2 ; z_dim++)
426 c_vector<int, DIM> neighbour;
427 neighbour(0) = x_dim;
428 neighbour(1) = y_dim;
429 neighbour(2) = z_dim;
431 neighbours.push_back(neighbour);
445 SetupLocalBoxes(neighbours);
454 for (
unsigned box_idx = 0 ; box_idx < mBoxes.size() ; box_idx++)
456 std::set<unsigned> local_boxes;
459 c_vector<unsigned, DIM> grid_indices = GetGridIndices(box_idx);
463 for (
unsigned neighbour = 0 ; neighbour < rNeighbours.size() ; neighbour++)
465 if (IsBoxInDomain(grid_indices + rNeighbours[neighbour]))
467 local_boxes.insert(GetLinearIndex(grid_indices + rNeighbours[neighbour]));
482 c_vector<bool, DIM> is_penultimate_periodic = IsIndexPenultimate(grid_indices);
483 unsigned num_penultimate_periodic_dimensions = 0;
485 for (
unsigned dim = 0 ; dim < DIM ; dim++)
487 is_penultimate_periodic(dim) = is_penultimate_periodic(dim) && mIsDomainPeriodic(dim);
488 if (is_penultimate_periodic(dim))
490 num_penultimate_periodic_dimensions++;
495 if (num_penultimate_periodic_dimensions > 0)
498 for (
unsigned dim = 0 ; dim < DIM ; dim++)
500 if (is_penultimate_periodic(dim))
505 c_vector<unsigned, DIM> ultimate_indices;
506 ultimate_indices = grid_indices;
507 ultimate_indices(dim) ++;
509 for (
unsigned neighbour = 0 ; neighbour < rNeighbours.size() ; neighbour++)
511 if (IsBoxInDomain(ultimate_indices + rNeighbours[neighbour]))
513 local_boxes.insert(GetLinearIndex(ultimate_indices + rNeighbours[neighbour]));
521 if (num_penultimate_periodic_dimensions > 1)
528 assert(mIsDomainPeriodic(0) && mIsDomainPeriodic(1));
530 local_boxes.insert(0);
539 if (is_penultimate_periodic(0) && is_penultimate_periodic(1))
541 c_vector<unsigned, DIM> ultimate_indices;
542 ultimate_indices = grid_indices;
543 ultimate_indices(0) ++;
544 ultimate_indices(1) ++;
546 for (
unsigned neighbour = 0 ; neighbour < rNeighbours.size() ; neighbour++)
548 if (IsBoxInDomain(ultimate_indices + rNeighbours[neighbour]))
550 local_boxes.insert(GetLinearIndex(ultimate_indices + rNeighbours[neighbour]));
556 if (is_penultimate_periodic(0) && is_penultimate_periodic(2))
558 c_vector<unsigned, DIM> ultimate_indices;
559 ultimate_indices = grid_indices;
560 ultimate_indices(0) ++;
561 ultimate_indices(2) ++;
563 for (
unsigned neighbour = 0 ; neighbour < rNeighbours.size() ; neighbour++)
565 if (IsBoxInDomain(ultimate_indices + rNeighbours[neighbour]))
567 local_boxes.insert(GetLinearIndex(ultimate_indices + rNeighbours[neighbour]));
573 if (is_penultimate_periodic(1) && is_penultimate_periodic(2))
575 c_vector<unsigned, DIM> ultimate_indices;
576 ultimate_indices = grid_indices;
577 ultimate_indices(1) ++;
578 ultimate_indices(2) ++;
580 for (
unsigned neighbour = 0 ; neighbour < rNeighbours.size() ; neighbour++)
582 if (IsBoxInDomain(ultimate_indices + rNeighbours[neighbour]))
584 local_boxes.insert(GetLinearIndex(ultimate_indices + rNeighbours[neighbour]));
590 if (num_penultimate_periodic_dimensions == 3)
592 assert(mIsDomainPeriodic(0) && mIsDomainPeriodic(1) && mIsDomainPeriodic(1));
593 local_boxes.insert(0);
606 mLocalBoxes.push_back(local_boxes);
633 for (
unsigned node_index = 0; node_index < rNodes.size(); node_index++)
635 rNodes[node_index]->ClearNeighbours();
637 unsigned box_index = CalculateContainingBox(rNodes[node_index]);
638 mBoxes[box_index].AddNode(rNodes[node_index]);
641 for (
unsigned i = 0; i < rNodes.size(); i++)
644 unsigned node_index = this_node->
GetIndex();
647 unsigned this_node_box_index = CalculateContainingBox(this_node);
650 std::set<unsigned> local_boxes_indices = rGetLocalBoxes(this_node_box_index);
653 for (std::set<unsigned>::iterator box_iter = local_boxes_indices.begin(); box_iter != local_boxes_indices.end();
657 std::set<Node<DIM>*>& r_contained_nodes = mBoxes[*box_iter].rGetNodesContained();
660 for (
typename std::set<
Node<DIM>*>::iterator node_iter = r_contained_nodes.begin();
661 node_iter != r_contained_nodes.end(); ++node_iter)
664 unsigned other_node_index = (*node_iter)->GetIndex();
667 if (*box_iter == this_node_box_index)
669 if (other_node_index > this_node->
GetIndex())
673 (*node_iter)->AddNeighbour(node_index);
680 (*node_iter)->AddNeighbour(node_index);
686 for (
unsigned i = 0; i < rNodes.size(); i++)
688 rNodes[i]->RemoveDuplicateNeighbours();