16 #ifndef dealii_cuda_hanging_nodes_internal_h 17 #define dealii_cuda_hanging_nodes_internal_h 19 #include <deal.II/base/config.h> 21 #ifdef DEAL_II_COMPILER_CUDA_AWARE 23 # include <deal.II/base/utilities.h> 25 # include <deal.II/dofs/dof_accessor.h> 26 # include <deal.II/dofs/dof_handler.h> 28 # include <deal.II/fe/fe_q.h> 29 # include <deal.II/fe/fe_tools.h> 31 DEAL_II_NAMESPACE_OPEN
54 const std::vector<unsigned int> &lexicographic_mapping);
59 template <
typename CellIterator>
62 std::vector<types::global_dof_index> & dof_indices,
63 const CellIterator & cell,
64 const std::shared_ptr<const Utilities::MPI::Partitioner> &partitioner,
65 unsigned int & mask)
const;
75 rotate_subface_index(
int times,
unsigned int &subface_index)
const;
78 rotate_face(
int times,
79 unsigned int n_dofs_1d,
80 std::vector<types::global_dof_index> &dofs)
const;
83 line_dof_idx(
int local_line,
85 unsigned int n_dofs_1d)
const;
88 transpose_face(std::vector<types::global_dof_index> &dofs)
const;
91 transpose_subface_index(
unsigned int &subface)
const;
94 using active_cell_iterator =
96 const unsigned int n_raw_lines;
97 std::vector<std::vector<std::pair<cell_iterator, unsigned int>>>
99 const std::vector<unsigned int> &lexicographic_mapping;
100 const unsigned int fe_degree;
107 # define DEAL_II_MAX_ELEM_DEGREE 10 108 __constant__
double constraint_weights[(DEAL_II_MAX_ELEM_DEGREE + 1) *
109 (DEAL_II_MAX_ELEM_DEGREE + 1)];
126 constexpr
unsigned int constr_type_x = 1 << 0;
127 constexpr
unsigned int constr_type_y = 1 << 1;
128 constexpr
unsigned int constr_type_z = 1 << 2;
131 constexpr
unsigned int constr_face_x = 1 << 3;
132 constexpr
unsigned int constr_face_y = 1 << 4;
133 constexpr
unsigned int constr_face_z = 1 << 5;
136 constexpr
unsigned int constr_edge_xy = 1 << 6;
137 constexpr
unsigned int constr_edge_yz = 1 << 7;
138 constexpr
unsigned int constr_edge_zx = 1 << 8;
142 setup_constraint_weigths(
unsigned int fe_degree)
147 fe_q.get_subface_interpolation_matrix(fe_q, 0, interpolation_matrix);
149 std::vector<unsigned int> mapping =
150 FETools::lexicographic_to_hierarchic_numbering<1>(
FE_Q<1>(fe_degree));
154 for (
unsigned int i = 0; i < fe_q.dofs_per_face; ++i)
155 for (
unsigned int j = 0; j < fe_q.dofs_per_face; ++j)
156 mapped_matrix(i, j) = interpolation_matrix(mapping[i], mapping[j]);
158 cudaError_t error_code =
159 cudaMemcpyToSymbol(internal::constraint_weights,
160 &mapped_matrix[0][0],
161 sizeof(
double) * fe_q.dofs_per_face *
171 unsigned int fe_degree,
173 const std::vector<unsigned int> &lexicographic_mapping)
174 : n_raw_lines(dof_handler.get_triangulation().n_raw_lines())
175 , line_to_cells(dim == 3 ? n_raw_lines : 0)
176 , lexicographic_mapping(lexicographic_mapping)
177 , fe_degree(fe_degree)
178 , dof_handler(dof_handler)
181 (dim == 3) || ((fe_degree % 2) == 1),
183 "This function is not implemented when dim = 2 and fe_degree is even."));
188 internal::setup_constraint_weigths<dim>(fe_degree);
210 const unsigned int line_to_children[12][2] = {{0, 2},
223 std::vector<std::vector<std::pair<cell_iterator, unsigned int>>>
224 line_to_inactive_cells(n_raw_lines);
227 for (
const auto &cell : dof_handler.cell_iterators())
229 for (
unsigned int line = 0; line < GeometryInfo<3>::lines_per_cell;
232 const unsigned int line_idx = cell->line(line)->index();
234 line_to_cells[line_idx].push_back(std::make_pair(cell, line));
236 line_to_inactive_cells[line_idx].push_back(
237 std::make_pair(cell, line));
244 for (
unsigned int line_idx = 0; line_idx < n_raw_lines; ++line_idx)
246 if ((line_to_cells[line_idx].size() > 0) &&
247 line_to_inactive_cells[line_idx].size() > 0)
251 const cell_iterator &inactive_cell =
252 line_to_inactive_cells[line_idx][0].first;
253 const unsigned int neighbor_line =
254 line_to_inactive_cells[line_idx][0].second;
256 for (
unsigned int c = 0; c < 2; ++c)
258 const cell_iterator &child =
259 inactive_cell->child(line_to_children[neighbor_line][c]);
260 const unsigned int child_line_idx =
261 child->line(neighbor_line)->index();
264 for (
const auto cl : line_to_cells[line_idx])
265 line_to_cells[child_line_idx].push_back(cl);
274 template <
typename CellIterator>
277 std::vector<types::global_dof_index> & dof_indices,
278 const CellIterator & cell,
279 const std::shared_ptr<const Utilities::MPI::Partitioner> &partitioner,
280 unsigned int & mask)
const 283 const unsigned int n_dofs_1d = fe_degree + 1;
284 const unsigned int dofs_per_face =
287 std::vector<types::global_dof_index> neighbor_dofs(dofs_per_face);
289 const auto lex_face_mapping =
291 FE_Q<dim - 1>(fe_degree));
293 for (
unsigned int face = 0; face < GeometryInfo<dim>::faces_per_cell;
296 if ((!cell->at_boundary(face)) &&
297 (cell->neighbor(face)->has_children() ==
false))
299 const active_cell_iterator &neighbor = cell->neighbor(face);
302 if (neighbor->level() < cell->level())
304 const unsigned int neighbor_face =
305 cell->neighbor_face_no(face);
308 unsigned int subface = 0;
309 for (; subface < GeometryInfo<dim>::max_children_per_face;
311 if (neighbor->neighbor_child_on_subface(neighbor_face,
316 neighbor->face(neighbor_face)->get_dof_indices(neighbor_dofs);
320 for (
auto &index : neighbor_dofs)
321 index = partitioner->global_to_local(index);
327 mask |= internal::constr_face_x;
329 mask |= internal::constr_type_x;
331 mask |= internal::constr_type_y;
335 mask |= internal::constr_face_y;
337 mask |= internal::constr_type_y;
339 mask |= internal::constr_type_x;
346 unsigned int offset = (face % 2 == 1) ? fe_degree : 0;
348 for (
unsigned int i = 0; i < n_dofs_1d; ++i)
350 unsigned int idx = 0;
353 idx = n_dofs_1d * offset + i;
356 idx = n_dofs_1d * i + offset;
358 dof_indices[idx] = neighbor_dofs[lex_face_mapping[i]];
363 const bool transpose = !(cell->face_orientation(face));
367 if (cell->face_rotation(face))
369 if (cell->face_flip(face))
372 rotate_face(rotate, n_dofs_1d, neighbor_dofs);
373 rotate_subface_index(rotate, subface);
377 transpose_face(neighbor_dofs);
378 transpose_subface_index(subface);
384 mask |= internal::constr_face_x;
386 mask |= internal::constr_type_x;
387 if (subface % 2 == 0)
388 mask |= internal::constr_type_y;
389 if (subface / 2 == 0)
390 mask |= internal::constr_type_z;
395 mask |= internal::constr_face_y;
397 mask |= internal::constr_type_y;
398 if (subface % 2 == 0)
399 mask |= internal::constr_type_z;
400 if (subface / 2 == 0)
401 mask |= internal::constr_type_x;
406 mask |= internal::constr_face_z;
408 mask |= internal::constr_type_z;
409 if (subface % 2 == 0)
410 mask |= internal::constr_type_x;
411 if (subface / 2 == 0)
412 mask |= internal::constr_type_y;
416 unsigned int offset = (face % 2 == 1) ? fe_degree : 0;
418 for (
unsigned int i = 0; i < n_dofs_1d; ++i)
420 for (
unsigned int j = 0; j < n_dofs_1d; ++j)
422 unsigned int idx = 0;
426 idx = n_dofs_1d * n_dofs_1d * i +
427 n_dofs_1d * j + offset;
431 idx = n_dofs_1d * n_dofs_1d * j +
432 n_dofs_1d * offset + i;
436 idx = n_dofs_1d * n_dofs_1d * offset +
440 neighbor_dofs[lex_face_mapping[n_dofs_1d * i +
458 const unsigned int line_to_edge[12][4] = {
459 {internal::constr_face_x | internal::constr_face_z,
460 internal::constr_edge_zx,
461 internal::constr_type_x | internal::constr_type_z,
462 internal::constr_type_y},
463 {internal::constr_face_x | internal::constr_face_z,
464 internal::constr_edge_zx,
465 internal::constr_type_z,
466 internal::constr_type_y},
467 {internal::constr_face_y | internal::constr_face_z,
468 internal::constr_edge_yz,
469 internal::constr_type_y | internal::constr_type_z,
470 internal::constr_type_x},
471 {internal::constr_face_y | internal::constr_face_z,
472 internal::constr_edge_yz,
473 internal::constr_type_z,
474 internal::constr_type_x},
475 {internal::constr_face_x | internal::constr_face_z,
476 internal::constr_edge_zx,
477 internal::constr_type_x,
478 internal::constr_type_y},
479 {internal::constr_face_x | internal::constr_face_z,
480 internal::constr_edge_zx,
482 internal::constr_type_y},
483 {internal::constr_face_y | internal::constr_face_z,
484 internal::constr_edge_yz,
485 internal::constr_type_y,
486 internal::constr_type_x},
487 {internal::constr_face_y | internal::constr_face_z,
488 internal::constr_edge_yz,
490 internal::constr_type_x},
491 {internal::constr_face_x | internal::constr_face_y,
492 internal::constr_edge_xy,
493 internal::constr_type_x | internal::constr_type_y,
494 internal::constr_type_z},
495 {internal::constr_face_x | internal::constr_face_y,
496 internal::constr_edge_xy,
497 internal::constr_type_y,
498 internal::constr_type_z},
499 {internal::constr_face_x | internal::constr_face_y,
500 internal::constr_edge_xy,
501 internal::constr_type_x,
502 internal::constr_type_z},
503 {internal::constr_face_x | internal::constr_face_y,
504 internal::constr_edge_xy,
506 internal::constr_type_z}};
508 for (
unsigned int local_line = 0;
509 local_line < GeometryInfo<dim>::lines_per_cell;
513 if (!(mask & line_to_edge[local_line][0]))
516 const unsigned int line = cell->line(local_line)->index();
517 for (
const auto edge_neighbor : line_to_cells[line])
520 const cell_iterator neighbor_cell = edge_neighbor.first;
521 if (neighbor_cell->level() < cell->level())
523 const unsigned int local_line_neighbor =
524 edge_neighbor.second;
525 mask |= line_to_edge[local_line][1] |
526 line_to_edge[local_line][2];
528 bool flipped =
false;
529 if (cell->line(local_line)->vertex_index(0) ==
530 neighbor_cell->line(local_line_neighbor)
535 mask |= line_to_edge[local_line][3];
537 else if (cell->line(local_line)->vertex_index(1) ==
538 neighbor_cell->line(local_line_neighbor)
543 else if (cell->line(local_line)->vertex_index(1) ==
544 neighbor_cell->line(local_line_neighbor)
550 else if (cell->line(local_line)->vertex_index(0) ==
551 neighbor_cell->line(local_line_neighbor)
555 mask |= line_to_edge[local_line][3];
562 neighbor_dofs.resize(n_dofs_1d * n_dofs_1d *
564 neighbor_cell->get_dof_indices(neighbor_dofs);
568 for (
auto &index : neighbor_dofs)
569 index = partitioner->global_to_local(index);
571 for (
unsigned int i = 0; i < n_dofs_1d; ++i)
574 const unsigned int idx =
575 line_dof_idx(local_line, i, n_dofs_1d);
576 dof_indices[idx] = neighbor_dofs
577 [lexicographic_mapping[line_dof_idx(
579 flipped ? fe_degree - i : i,
597 unsigned int &subface_index)
const 599 const unsigned int rot_mapping[4] = {2, 0, 3, 1};
602 times = times < 0 ? times + 4 : times;
603 for (
int t = 0; t < times; ++t)
604 subface_index = rot_mapping[subface_index];
611 HangingNodes<dim>::rotate_face(
613 unsigned int n_dofs_1d,
614 std::vector<types::global_dof_index> &dofs)
const 616 const unsigned int rot_mapping[4] = {2, 0, 3, 1};
619 times = times < 0 ? times + 4 : times;
621 std::vector<types::global_dof_index> copy(dofs.size());
622 for (
int t = 0; t < times; ++t)
627 for (
unsigned int i = 0; i < 4; ++i)
628 dofs[rot_mapping[i]] = copy[i];
631 const unsigned int n_int = n_dofs_1d - 2;
632 unsigned int offset = 4;
633 for (
unsigned int i = 0; i < n_int; ++i)
636 dofs[offset + i] = copy[offset + 2 * n_int + (n_int - 1 - i)];
638 dofs[offset + n_int + i] =
639 copy[offset + 3 * n_int + (n_int - 1 - i)];
641 dofs[offset + 2 * n_int + i] = copy[offset + n_int + i];
643 dofs[offset + 3 * n_int + i] = copy[offset + i];
649 for (
unsigned int i = 0; i < n_int; ++i)
650 for (
unsigned int j = 0; j < n_int; ++j)
651 dofs[offset + i * n_int + j] =
652 copy[offset + j * n_int + (n_int - 1 - i)];
660 HangingNodes<dim>::line_dof_idx(
int local_line,
662 unsigned int n_dofs_1d)
const 664 unsigned int x, y, z;
669 (local_line % 4 == 0) ? 0 : (local_line % 4 == 1) ? fe_degree : dof;
671 (local_line % 4 == 2) ? 0 : (local_line % 4 == 3) ? fe_degree : dof;
672 z = (local_line / 4) * fe_degree;
676 x = ((local_line - 8) % 2) * fe_degree;
677 y = ((local_line - 8) / 2) * fe_degree;
681 return n_dofs_1d * n_dofs_1d * z + n_dofs_1d * y + x;
688 HangingNodes<dim>::transpose_face(
689 std::vector<types::global_dof_index> &dofs)
const 691 const std::vector<types::global_dof_index> copy(dofs);
698 const unsigned int n_int = fe_degree - 1;
699 unsigned int offset = 4;
700 for (
unsigned int i = 0; i < n_int; ++i)
703 dofs[offset + i] = copy[offset + 2 * n_int + i];
705 dofs[offset + n_int + i] = copy[offset + 3 * n_int + i];
707 dofs[offset + 2 * n_int + i] = copy[offset + i];
709 dofs[offset + 3 * n_int + i] = copy[offset + n_int + i];
714 for (
unsigned int i = 0; i < n_int; ++i)
715 for (
unsigned int j = 0; j < n_int; ++j)
716 dofs[offset + i * n_int + j] = copy[offset + j * n_int + i];
723 HangingNodes<dim>::transpose_subface_index(
unsigned int &subface)
const 727 else if (subface == 2)
738 template <
unsigned int size>
739 __device__
inline unsigned int 740 index2(
unsigned int i,
unsigned int j)
747 template <
unsigned int size>
748 __device__
inline unsigned int 749 index3(
unsigned int i,
unsigned int j,
unsigned int k)
751 return i + size * j + size * size * k;
756 template <
unsigned int fe_degree,
757 unsigned int direction,
760 __device__
inline void 761 interpolate_boundary_2d(
const unsigned int constr, Number *values)
763 const unsigned int x_idx = threadIdx.x % (fe_degree + 1);
764 const unsigned int y_idx = threadIdx.y;
766 const unsigned int this_type =
767 (direction == 0) ? internal::constr_type_x : internal::constr_type_y;
769 if (constr & (((direction == 0) ? internal::constr_face_y : 0) |
770 ((direction == 1) ? internal::constr_face_x : 0)))
772 const unsigned int interp_idx = (direction == 0) ? x_idx : y_idx;
779 ((constr & internal::constr_type_y) ? (y_idx == 0) :
780 (y_idx == fe_degree))) ||
782 ((constr & internal::constr_type_x) ? (x_idx == 0) :
783 (x_idx == fe_degree)));
787 const bool type = constr & this_type;
791 for (
unsigned int i = 0; i <= fe_degree; ++i)
793 const unsigned int real_idx =
794 (direction == 0) ? index2<fe_degree + 1>(i, y_idx) :
795 index2<fe_degree + 1>(x_idx, i);
799 internal::constraint_weights[i * (fe_degree + 1) +
801 internal::constraint_weights[interp_idx *
804 t +=
w * values[real_idx];
809 for (
unsigned int i = 0; i <= fe_degree; ++i)
811 const unsigned int real_idx =
812 (direction == 0) ? index2<fe_degree + 1>(i, y_idx) :
813 index2<fe_degree + 1>(x_idx, i);
817 internal::constraint_weights[(fe_degree - i) *
821 internal::constraint_weights[(fe_degree -
825 t +=
w * values[real_idx];
833 values[index2<fe_degree + 1>(x_idx, y_idx)] = t;
839 template <
unsigned int fe_degree,
840 unsigned int direction,
843 __device__
inline void 844 interpolate_boundary_3d(
const unsigned int constr, Number *values)
846 const unsigned int x_idx = threadIdx.x % (fe_degree + 1);
847 const unsigned int y_idx = threadIdx.y;
848 const unsigned int z_idx = threadIdx.z;
850 const unsigned int this_type =
851 (direction == 0) ? internal::constr_type_x :
852 (direction == 1) ? internal::constr_type_y :
853 internal::constr_type_z;
854 const unsigned int face1_type =
855 (direction == 0) ? internal::constr_type_y :
856 (direction == 1) ? internal::constr_type_z :
857 internal::constr_type_x;
858 const unsigned int face2_type =
859 (direction == 0) ? internal::constr_type_z :
860 (direction == 1) ? internal::constr_type_x :
861 internal::constr_type_y;
865 const unsigned int face1 = (direction == 0) ? internal::constr_face_y :
867 internal::constr_face_z :
868 internal::constr_face_x;
869 const unsigned int face2 = (direction == 0) ? internal::constr_face_z :
871 internal::constr_face_x :
872 internal::constr_face_y;
873 const unsigned int edge = (direction == 0) ? internal::constr_edge_yz :
875 internal::constr_edge_zx :
876 internal::constr_edge_xy;
878 if (constr & (face1 | face2 | edge))
880 const unsigned int interp_idx =
881 (direction == 0) ? x_idx : (direction == 1) ? y_idx : z_idx;
882 const unsigned int face1_idx =
883 (direction == 0) ? y_idx : (direction == 1) ? z_idx : x_idx;
884 const unsigned int face2_idx =
885 (direction == 0) ? z_idx : (direction == 1) ? x_idx : y_idx;
890 const bool on_face1 = (constr & face1_type) ?
892 (face1_idx == fe_degree);
893 const bool on_face2 = (constr & face2_type) ?
895 (face2_idx == fe_degree);
896 const bool flag = (((constr & face1) && on_face1) ||
897 ((constr & face2) && on_face2) ||
898 ((constr & edge) && on_face1 && on_face2));
902 const bool type = constr & this_type;
905 for (
unsigned int i = 0; i <= fe_degree; ++i)
907 const unsigned int real_idx =
909 index3<fe_degree + 1>(i, y_idx, z_idx) :
911 index3<fe_degree + 1>(x_idx, i, z_idx) :
912 index3<fe_degree + 1>(x_idx, y_idx, i);
916 internal::constraint_weights[i * (fe_degree + 1) +
918 internal::constraint_weights[interp_idx *
921 t +=
w * values[real_idx];
926 for (
unsigned int i = 0; i <= fe_degree; ++i)
928 const unsigned int real_idx =
930 index3<fe_degree + 1>(i, y_idx, z_idx) :
932 index3<fe_degree + 1>(x_idx, i, z_idx) :
933 index3<fe_degree + 1>(x_idx, y_idx, i);
937 internal::constraint_weights[(fe_degree - i) *
941 internal::constraint_weights[(fe_degree -
945 t +=
w * values[real_idx];
953 values[index3<fe_degree + 1>(x_idx, y_idx, z_idx)] = t;
968 template <
int dim,
int fe_degree,
bool transpose,
typename Number>
970 resolve_hanging_nodes(
const unsigned int constr, Number *values)
974 internal::interpolate_boundary_2d<fe_degree, 0, transpose>(constr,
976 internal::interpolate_boundary_2d<fe_degree, 1, transpose>(constr,
982 internal::interpolate_boundary_3d<fe_degree, 0, transpose>(constr,
985 internal::interpolate_boundary_3d<fe_degree, 1, transpose>(constr,
988 internal::interpolate_boundary_3d<fe_degree, 2, transpose>(constr,
995 DEAL_II_NAMESPACE_CLOSE
void setup_line_to_cell()
#define AssertThrow(cond, exc)
static ::ExceptionBase & ExcMessage(std::string arg1)
HangingNodes(unsigned int fe_degree, const DoFHandler< dim > &dof_handler, const std::vector< unsigned int > &lexicographic_mapping)
#define AssertCuda(error_code)
void setup_constraints(std::vector< types::global_dof_index > &dof_indices, const CellIterator &cell, const std::shared_ptr< const Utilities::MPI::Partitioner > &partitioner, unsigned int &mask) const
SymmetricTensor< rank_, dim, Number > transpose(const SymmetricTensor< rank_, dim, Number > &t)
Tensor< 2, dim, Number > w(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
void swap(Vector< Number > &u, Vector< Number > &v)
static ::ExceptionBase & ExcNotImplemented()
static ::ExceptionBase & ExcInternalError()