17 #include <deal.II/distributed/tria.h> 19 #include <deal.II/dofs/dof_tools.h> 21 #include <deal.II/fe/fe_tools.h> 23 #include <deal.II/matrix_free/shape_info.h> 25 #include <deal.II/multigrid/mg_transfer_internal.h> 27 DEAL_II_NAMESPACE_OPEN
42 DoFPair(
const unsigned int level,
46 , global_dof_index(global_dof_index)
47 , level_dof_index(level_dof_index)
51 : level(
numbers::invalid_unsigned_int)
59 template <
int dim,
int spacedim>
62 const ::DoFHandler<dim, spacedim> &mg_dof,
64 std::vector<std::vector<
65 std::pair<types::global_dof_index, types::global_dof_index>>>
67 std::vector<std::vector<
68 std::pair<types::global_dof_index, types::global_dof_index>>>
69 ©_indices_global_mine,
70 std::vector<std::vector<
71 std::pair<types::global_dof_index, types::global_dof_index>>>
72 & copy_indices_level_mine,
73 const bool skip_interface_dofs)
86 std::vector<DoFPair> send_data_temp;
88 const unsigned int n_levels =
89 mg_dof.get_triangulation().n_global_levels();
90 copy_indices.resize(n_levels);
91 copy_indices_global_mine.resize(n_levels);
92 copy_indices_level_mine.resize(n_levels);
96 const unsigned int dofs_per_cell = mg_dof.get_fe().dofs_per_cell;
97 std::vector<types::global_dof_index> global_dof_indices(dofs_per_cell);
98 std::vector<types::global_dof_index> level_dof_indices(dofs_per_cell);
100 for (
unsigned int level = 0; level < n_levels; ++level)
102 std::vector<bool> dof_touched(globally_relevant.
n_elements(),
false);
103 copy_indices[level].clear();
104 copy_indices_level_mine[level].clear();
105 copy_indices_global_mine[level].clear();
107 typename ::DoFHandler<dim, spacedim>::active_cell_iterator
108 level_cell = mg_dof.begin_active(level);
109 const typename ::DoFHandler<dim, spacedim>::active_cell_iterator
110 level_end = mg_dof.end_active(level);
112 for (; level_cell != level_end; ++level_cell)
114 if (mg_dof.get_triangulation().locally_owned_subdomain() !=
116 (level_cell->level_subdomain_id() ==
118 level_cell->subdomain_id() ==
124 level_cell->get_dof_indices(global_dof_indices);
125 level_cell->get_mg_dof_indices(level_dof_indices);
127 for (
unsigned int i = 0; i < dofs_per_cell; ++i)
131 if (skip_interface_dofs && mg_constrained_dofs !=
nullptr &&
133 level, level_dof_indices[i]))
140 if (dof_touched[global_idx])
142 bool global_mine = mg_dof.locally_owned_dofs().is_element(
143 global_dof_indices[i]);
145 mg_dof.locally_owned_mg_dofs(level).is_element(
146 level_dof_indices[i]);
149 if (global_mine && level_mine)
151 copy_indices[level].emplace_back(global_dof_indices[i],
152 level_dof_indices[i]);
154 else if (global_mine)
156 copy_indices_global_mine[level].emplace_back(
157 global_dof_indices[i], level_dof_indices[i]);
160 send_data_temp.emplace_back(level,
161 global_dof_indices[i],
162 level_dof_indices[i]);
169 dof_touched[global_idx] =
true;
174 const ::parallel::Triangulation<dim, spacedim> *tria =
176 &mg_dof.get_triangulation()));
178 send_data_temp.size() == 0 || tria !=
nullptr,
180 "We should only be sending information with a parallel Triangulation!"));
182 #ifdef DEAL_II_WITH_MPI 190 const std::set<types::subdomain_id> &neighbors =
191 tria->level_ghost_owners();
192 std::map<int, std::vector<DoFPair>> send_data;
196 for (
typename std::vector<DoFPair>::iterator dofpair =
197 send_data_temp.begin();
198 dofpair != send_data_temp.end();
201 std::set<types::subdomain_id>::iterator it;
202 for (it = neighbors.begin(); it != neighbors.end(); ++it)
205 .locally_owned_mg_dofs_per_processor(
207 .is_element(dofpair->level_dof_index))
209 send_data[*it].push_back(*dofpair);
215 Assert(it != neighbors.end(),
220 std::vector<MPI_Request> requests;
222 for (std::set<types::subdomain_id>::iterator it = neighbors.begin();
223 it != neighbors.end();
226 requests.push_back(MPI_Request());
227 unsigned int dest = *it;
228 std::vector<DoFPair> &data = send_data[dest];
234 const int ierr = MPI_Isend(data.data(),
235 data.size() *
sizeof(data[0]),
239 tria->get_communicator(),
240 &*requests.rbegin());
245 const int ierr = MPI_Isend(
nullptr,
250 tria->get_communicator(),
251 &*requests.rbegin());
260 std::vector<DoFPair> receive_buffer;
261 for (
unsigned int counter = 0; counter < neighbors.size();
266 int ierr = MPI_Probe(MPI_ANY_SOURCE,
268 tria->get_communicator(),
271 ierr = MPI_Get_count(&status, MPI_BYTE, &len);
276 ierr = MPI_Recv(
nullptr,
281 tria->get_communicator(),
287 int count = len /
sizeof(DoFPair);
288 Assert(static_cast<int>(count *
sizeof(DoFPair)) == len,
290 receive_buffer.resize(count);
292 void *ptr = receive_buffer.data();
298 tria->get_communicator(),
302 for (
unsigned int i = 0; i < receive_buffer.size(); ++i)
304 copy_indices_level_mine[receive_buffer[i].level]
305 .emplace_back(receive_buffer[i].global_dof_index,
306 receive_buffer[i].level_dof_index);
312 if (requests.size() > 0)
314 const int ierr = MPI_Waitall(requests.size(),
316 MPI_STATUSES_IGNORE);
324 const int ierr = MPI_Barrier(tria->get_communicator());
333 std::less<std::pair<types::global_dof_index, types::global_dof_index>>
335 for (
unsigned int level = 0; level < copy_indices.size(); ++level)
336 std::sort(copy_indices[level].begin(),
337 copy_indices[level].end(),
339 for (
unsigned int level = 0; level < copy_indices_level_mine.size();
341 std::sort(copy_indices_level_mine[level].begin(),
342 copy_indices_level_mine[level].end(),
344 for (
unsigned int level = 0; level < copy_indices_global_mine.size();
346 std::sort(copy_indices_global_mine[level].begin(),
347 copy_indices_global_mine[level].end(),
355 template <
typename Number>
357 reinit_ghosted_vector(
359 std::vector<types::global_dof_index> & ghosted_level_dofs,
360 const MPI_Comm & communicator,
362 std::vector<std::pair<unsigned int, unsigned int>>
363 ©_indices_global_mine)
365 std::sort(ghosted_level_dofs.begin(), ghosted_level_dofs.end());
367 ghosted_dofs.
add_indices(ghosted_level_dofs.begin(),
368 std::unique(ghosted_level_dofs.begin(),
369 ghosted_level_dofs.end()));
370 ghosted_dofs.compress();
373 if (ghosted_level_vector.
size() == locally_owned.
size())
378 ghosted_dofs.add_indices(part->ghost_indices());
379 for (
unsigned int i = 0; i < copy_indices_global_mine.size(); ++i)
380 copy_indices_global_mine[i].second =
382 ghosted_dofs.index_within_set(
383 part->local_to_global(copy_indices_global_mine[i].second));
385 ghosted_level_vector.
reinit(locally_owned, ghosted_dofs, communicator);
390 copy_indices_to_mpi_local_numbers(
392 const std::vector<types::global_dof_index> &mine,
393 const std::vector<types::global_dof_index> &remote,
394 std::vector<unsigned int> & localized_indices)
396 localized_indices.resize(mine.size() + remote.size(),
398 for (
unsigned int i = 0; i < mine.size(); ++i)
402 for (
unsigned int i = 0; i < remote.size(); ++i)
411 compute_shift_within_children(
const unsigned int child,
412 const unsigned int fe_shift_1d,
413 const unsigned int fe_degree)
417 unsigned int c_tensor_index[dim];
418 unsigned int tmp = child;
419 for (
unsigned int d = 0;
d < dim; ++
d)
421 c_tensor_index[
d] = tmp % 2;
424 const unsigned int n_child_dofs_1d = fe_degree + 1 + fe_shift_1d;
425 unsigned int factor = 1;
426 unsigned int shift = fe_shift_1d * c_tensor_index[0];
427 for (
unsigned int d = 1;
d < dim; ++
d)
429 factor *= n_child_dofs_1d;
430 shift = shift + factor * fe_shift_1d * c_tensor_index[
d];
440 const unsigned int child,
441 const unsigned int fe_shift_1d,
442 const unsigned int fe_degree,
443 const std::vector<unsigned int> & lexicographic_numbering,
444 const std::vector<types::global_dof_index> &local_dof_indices,
447 const unsigned int n_child_dofs_1d = fe_degree + 1 + fe_shift_1d;
448 const unsigned int shift =
449 compute_shift_within_children<dim>(child, fe_shift_1d, fe_degree);
451 local_dof_indices.size() / Utilities::fixed_power<dim>(fe_degree + 1);
453 const unsigned int n_scalar_cell_dofs =
454 Utilities::fixed_power<dim>(n_child_dofs_1d);
455 for (
unsigned int c = 0, m = 0; c < n_components; ++c)
456 for (
unsigned int k = 0; k < (dim > 2 ? (fe_degree + 1) : 1); ++k)
457 for (
unsigned int j = 0; j < (dim > 1 ? (fe_degree + 1) : 1); ++j)
458 for (
unsigned int i = 0; i < (fe_degree + 1); ++i, ++m)
460 const unsigned int index =
461 c * n_scalar_cell_dofs +
462 k * n_child_dofs_1d * n_child_dofs_1d + j * n_child_dofs_1d +
466 local_dof_indices[lexicographic_numbering[m]],
468 indices[index] = local_dof_indices[lexicographic_numbering[m]];
472 template <
int dim,
typename Number>
474 setup_element_info(ElementInfo<Number> & elem_info,
476 const ::DoFHandler<dim> &mg_dof)
479 elem_info.n_components = mg_dof.get_fe().element_multiplicity(0);
481 elem_info.n_components,
482 mg_dof.get_fe().dofs_per_cell);
484 elem_info.fe_degree = fe.
degree;
509 const unsigned int n_child_dofs_1d =
513 elem_info.n_child_cell_dofs =
514 elem_info.n_components * Utilities::fixed_power<dim>(n_child_dofs_1d);
518 shape_info.
reinit(dummy_quadrature, mg_dof.get_fe(), 0);
525 for (
unsigned int c = 0; c < GeometryInfo<1>::max_children_per_cell; ++c)
529 .prolongation_matrix_1d[i * n_child_dofs_1d + j + c * shift] =
542 const unsigned int level,
543 std::vector<types::global_dof_index> &dof_indices)
545 if (mg_constrained_dofs !=
nullptr &&
548 for (
auto &ind : dof_indices)
566 template <
int dim,
typename Number>
569 const ::DoFHandler<dim> & mg_dof,
571 ElementInfo<Number> & elem_info,
572 std::vector<std::vector<unsigned int>> &level_dof_indices,
573 std::vector<std::vector<std::pair<unsigned int, unsigned int>>>
574 & parent_child_connect,
575 std::vector<unsigned int> &n_owned_level_cells,
576 std::vector<std::vector<std::vector<unsigned short>>> &dirichlet_indices,
577 std::vector<std::vector<Number>> & weights_on_refined,
578 std::vector<std::vector<std::pair<unsigned int, unsigned int>>>
579 ©_indices_global_mine,
581 &ghosted_level_vector)
583 level_dof_indices.clear();
584 parent_child_connect.clear();
585 n_owned_level_cells.clear();
586 dirichlet_indices.clear();
587 weights_on_refined.clear();
593 const ::Triangulation<dim> &tria = mg_dof.get_triangulation();
599 std::string fe_name = mg_dof.get_fe().base_element(0).get_name();
601 const std::size_t template_starts = fe_name.find_first_of(
'<');
602 Assert(fe_name[template_starts + 1] ==
603 (dim == 1 ?
'1' : (dim == 2 ?
'2' :
'3')),
605 fe_name[template_starts + 1] =
'1';
607 const std::unique_ptr<FiniteElement<1>> fe(
608 FETools::get_fe_by_name<1, 1>(fe_name));
610 setup_element_info(elem_info, *fe, mg_dof);
615 const unsigned int n_levels = tria.n_global_levels();
616 level_dof_indices.resize(n_levels);
617 parent_child_connect.resize(n_levels - 1);
618 n_owned_level_cells.resize(n_levels - 1);
619 std::vector<std::vector<unsigned int>> coarse_level_indices(n_levels - 1);
620 for (
unsigned int level = 0;
621 level < std::min(tria.n_levels(), n_levels - 1);
623 coarse_level_indices[level].resize(tria.n_raw_cells(level),
625 std::vector<types::global_dof_index> local_dof_indices(
626 mg_dof.get_fe().dofs_per_cell);
627 dirichlet_indices.resize(n_levels - 1);
632 if (ghosted_level_vector.max_level() != n_levels - 1)
633 ghosted_level_vector.resize(0, n_levels - 1);
635 for (
unsigned int level = n_levels - 1; level > 0; --level)
637 unsigned int counter = 0;
638 std::vector<types::global_dof_index> global_level_dof_indices;
639 std::vector<types::global_dof_index> global_level_dof_indices_remote;
640 std::vector<types::global_dof_index> ghosted_level_dofs;
641 std::vector<types::global_dof_index> global_level_dof_indices_l0;
642 std::vector<types::global_dof_index> ghosted_level_dofs_l0;
645 typename ::DoFHandler<dim>::cell_iterator cell,
646 endc = mg_dof.end(level - 1);
647 for (cell = mg_dof.begin(level - 1); cell != endc; ++cell)
651 if (!cell->has_children())
654 bool consider_cell =
false;
655 if (tria.locally_owned_subdomain() ==
657 cell->level_subdomain_id() == tria.locally_owned_subdomain())
658 consider_cell =
true;
663 bool cell_is_remote = !consider_cell;
664 for (
unsigned int c = 0;
665 c < GeometryInfo<dim>::max_children_per_cell;
667 if (cell->child(c)->level_subdomain_id() ==
668 tria.locally_owned_subdomain())
670 consider_cell =
true;
685 std::vector<types::global_dof_index> &next_indices =
686 cell_is_remote ? global_level_dof_indices_remote :
687 global_level_dof_indices;
688 const std::size_t start_index = next_indices.size();
689 next_indices.resize(start_index + elem_info.n_child_cell_dofs,
691 for (
unsigned int c = 0;
692 c < GeometryInfo<dim>::max_children_per_cell;
695 if (cell_is_remote && cell->child(c)->level_subdomain_id() !=
696 tria.locally_owned_subdomain())
698 cell->child(c)->get_mg_dof_indices(local_dof_indices);
700 replace(mg_constrained_dofs, level, local_dof_indices);
703 mg_dof.locally_owned_mg_dofs(level);
704 for (
unsigned int i = 0; i < local_dof_indices.size(); ++i)
705 if (!owned_level_dofs.
is_element(local_dof_indices[i]))
706 ghosted_level_dofs.push_back(local_dof_indices[i]);
708 add_child_indices<dim>(c,
712 elem_info.lexicographic_numbering,
714 &next_indices[start_index]);
717 if (cell->child(c)->has_children() &&
718 (tria.locally_owned_subdomain() ==
720 cell->child(c)->level_subdomain_id() ==
721 tria.locally_owned_subdomain()))
723 const unsigned int child_index =
724 coarse_level_indices[level][cell->child(c)->index()];
726 parent_child_connect[level].size());
727 unsigned int parent_index = counter;
736 start_index / elem_info.n_child_cell_dofs +
738 parent_child_connect[level][child_index] =
739 std::make_pair(parent_index, c);
741 static_cast<unsigned short>(-1));
745 if (mg_constrained_dofs !=
nullptr)
746 for (
unsigned int i = 0;
747 i < mg_dof.get_fe().dofs_per_cell;
752 [elem_info.lexicographic_numbering[i]]))
753 dirichlet_indices[level][child_index].push_back(i);
759 coarse_level_indices[level - 1].size());
760 coarse_level_indices[level - 1][cell->index()] = counter++;
767 if (level == 1 && !cell_is_remote)
769 cell->get_mg_dof_indices(local_dof_indices);
771 replace(mg_constrained_dofs, level - 1, local_dof_indices);
773 const IndexSet &owned_level_dofs_l0 =
774 mg_dof.locally_owned_mg_dofs(0);
775 for (
unsigned int i = 0; i < local_dof_indices.size(); ++i)
776 if (!owned_level_dofs_l0.
is_element(local_dof_indices[i]))
777 ghosted_level_dofs_l0.push_back(local_dof_indices[i]);
779 const std::size_t start_index =
780 global_level_dof_indices_l0.size();
781 global_level_dof_indices_l0.resize(
782 start_index + elem_info.n_child_cell_dofs,
784 add_child_indices<dim>(
788 elem_info.lexicographic_numbering,
790 &global_level_dof_indices_l0[start_index]);
792 dirichlet_indices[0].emplace_back();
793 if (mg_constrained_dofs !=
nullptr)
794 for (
unsigned int i = 0; i < mg_dof.get_fe().dofs_per_cell;
798 local_dof_indices[elem_info
799 .lexicographic_numbering[i]]))
800 dirichlet_indices[0].back().push_back(i);
808 global_level_dof_indices.size());
809 n_owned_level_cells[level - 1] = counter;
810 dirichlet_indices[level - 1].resize(counter);
811 parent_child_connect[level - 1].resize(
819 if (level < n_levels - 1)
820 for (std::vector<std::pair<unsigned int, unsigned int>>::iterator
821 i = parent_child_connect[level].begin();
822 i != parent_child_connect[level].end();
824 if (i->first >= tria.n_cells(level))
826 i->first -= tria.n_cells(level);
833 const MPI_Comm communicator =
836 reinit_ghosted_vector(mg_dof.locally_owned_mg_dofs(level),
839 ghosted_level_vector[level],
840 copy_indices_global_mine[level]);
842 copy_indices_to_mpi_local_numbers(
843 *ghosted_level_vector[level].get_partitioner(),
844 global_level_dof_indices,
845 global_level_dof_indices_remote,
846 level_dof_indices[level]);
850 for (
unsigned int i = 0; i < parent_child_connect[0].size(); ++i)
851 parent_child_connect[0][i] = std::make_pair(i, 0U);
853 reinit_ghosted_vector(mg_dof.locally_owned_mg_dofs(0),
854 ghosted_level_dofs_l0,
856 ghosted_level_vector[0],
857 copy_indices_global_mine[0]);
859 copy_indices_to_mpi_local_numbers(
860 *ghosted_level_vector[0].get_partitioner(),
861 global_level_dof_indices_l0,
862 std::vector<types::global_dof_index>(),
863 level_dof_indices[0]);
869 const unsigned int n_child_dofs_1d =
874 weights_on_refined.resize(n_levels - 1);
875 for (
unsigned int level = 1; level < n_levels; ++level)
877 ghosted_level_vector[level] = 0;
878 for (
unsigned int c = 0; c < n_owned_level_cells[level - 1]; ++c)
879 for (
unsigned int j = 0; j < elem_info.n_child_cell_dofs; ++j)
880 ghosted_level_vector[level].local_element(
881 level_dof_indices[level][elem_info.n_child_cell_dofs * c +
884 ghosted_level_vector[level].update_ghost_values();
886 std::vector<unsigned int> degree_to_3(n_child_dofs_1d);
888 for (
unsigned int i = 1; i < n_child_dofs_1d - 1; ++i)
890 degree_to_3.back() = 2;
894 weights_on_refined[level - 1].resize(n_owned_level_cells[level - 1] *
895 Utilities::fixed_power<dim>(3));
896 for (
unsigned int c = 0; c < n_owned_level_cells[level - 1]; ++c)
897 for (
unsigned int k = 0, m = 0; k < (dim > 2 ? n_child_dofs_1d : 1);
899 for (
unsigned int j = 0; j < (dim > 1 ? n_child_dofs_1d : 1); ++j)
901 unsigned int shift = 9 * degree_to_3[k] + 3 * degree_to_3[j];
902 for (
unsigned int i = 0; i < n_child_dofs_1d; ++i, ++m)
903 weights_on_refined[level -
904 1][c * Utilities::fixed_power<dim>(3) +
905 shift + degree_to_3[i]] =
907 ghosted_level_vector[level].local_element(
908 level_dof_indices[level]
909 [elem_info.n_child_cell_dofs * c + m]);
919 #include "mg_transfer_internal.inst" 921 DEAL_II_NAMESPACE_CLOSE
bool is_boundary_index(const unsigned int level, const types::global_dof_index index) const
void reinit(const size_type size, const bool omit_zeroing_entries=false)
const std::shared_ptr< const Utilities::MPI::Partitioner > & get_partitioner() const
static const unsigned int invalid_unsigned_int
const std::vector< Point< dim > > & get_unit_support_points() const
bool at_refinement_edge(const unsigned int level, const types::global_dof_index index) const
const AffineConstraints< double > & get_level_constraint_matrix(const unsigned int level) const
#define AssertDimension(dim1, dim2)
const types::subdomain_id invalid_subdomain_id
void reinit(const Quadrature< 1 > &quad, const FiniteElement< dim > &fe_dim, const unsigned int base_element=0)
#define AssertIndexRange(index, range)
size_type n_elements() const
void add_indices(const ForwardIterator &begin, const ForwardIterator &end)
const unsigned int degree
#define AssertThrow(cond, exc)
const unsigned int dofs_per_line
unsigned long long int global_dof_index
virtual MPI_Comm get_communicator() const
virtual const FullMatrix< double > & get_prolongation_matrix(const unsigned int child, const RefinementCase< dim > &refinement_case=RefinementCase< dim >::isotropic_refinement) const
static::ExceptionBase & ExcMessage(std::string arg1)
virtual size_type size() const override
#define Assert(cond, exc)
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
unsigned int global_to_local(const types::global_dof_index global_index) const
const unsigned int dofs_per_cell
const types::subdomain_id artificial_subdomain_id
#define AssertThrowMPI(error_code)
bool is_identity_constrained(const size_type index) const
size_type index_within_set(const size_type global_index) const
bool is_element(const size_type index) const
static::ExceptionBase & ExcNotImplemented()
const unsigned int dofs_per_vertex
const std::vector< std::pair< size_type, number > > * get_constraint_entries(const size_type line) const
std::vector< unsigned int > lexicographic_numbering
const types::global_dof_index invalid_dof_index
static::ExceptionBase & ExcInternalError()
size_type n_constraints() const