17 #include <deal.II/distributed/tria.h> 18 #include <deal.II/dofs/dof_tools.h> 19 #include <deal.II/fe/fe_tools.h> 20 #include <deal.II/matrix_free/shape_info.h> 21 #include <deal.II/multigrid/mg_transfer_internal.h> 23 DEAL_II_NAMESPACE_OPEN
39 DoFPair(
const unsigned int level,
45 level_dof_index(level_dof_index)
50 level (
numbers::invalid_unsigned_int),
59 template <
int dim,
int spacedim>
60 void fill_copy_indices(const ::DoFHandler<dim,spacedim> &mg_dof,
62 std::vector<std::vector<std::pair<types::global_dof_index, types::global_dof_index> > > ©_indices,
63 std::vector<std::vector<std::pair<types::global_dof_index, types::global_dof_index> > > ©_indices_global_mine,
64 std::vector<std::vector<std::pair<types::global_dof_index, types::global_dof_index> > > ©_indices_level_mine,
65 const bool skip_interface_dofs)
78 std::vector<DoFPair> send_data_temp;
80 const unsigned int n_levels = mg_dof.get_triangulation().n_global_levels();
81 copy_indices.resize(n_levels);
82 copy_indices_global_mine.resize(n_levels);
83 copy_indices_level_mine.resize(n_levels);
87 const unsigned int dofs_per_cell = mg_dof.get_fe().dofs_per_cell;
88 std::vector<types::global_dof_index> global_dof_indices (dofs_per_cell);
89 std::vector<types::global_dof_index> level_dof_indices (dofs_per_cell);
91 for (
unsigned int level=0; level<n_levels; ++level)
93 std::vector<bool> dof_touched(globally_relevant.
n_elements(),
false);
94 copy_indices[level].clear();
95 copy_indices_level_mine[level].clear();
96 copy_indices_global_mine[level].clear();
98 typename ::DoFHandler<dim,spacedim>::active_cell_iterator
99 level_cell = mg_dof.begin_active(level);
100 const typename ::DoFHandler<dim,spacedim>::active_cell_iterator
101 level_end = mg_dof.end_active(level);
103 for (; level_cell!=level_end; ++level_cell)
113 level_cell->get_dof_indices (global_dof_indices);
114 level_cell->get_mg_dof_indices (level_dof_indices);
116 for (
unsigned int i=0; i<dofs_per_cell; ++i)
119 if (skip_interface_dofs &&
120 mg_constrained_dofs !=
nullptr 126 if (dof_touched[global_idx])
128 bool global_mine = mg_dof.locally_owned_dofs().is_element(global_dof_indices[i]);
129 bool level_mine = mg_dof.locally_owned_mg_dofs(level).is_element(level_dof_indices[i]);
132 if (global_mine && level_mine)
134 copy_indices[level].emplace_back(global_dof_indices[i],
135 level_dof_indices[i]);
137 else if (global_mine)
139 copy_indices_global_mine[level].emplace_back(global_dof_indices[i],
140 level_dof_indices[i]);
143 send_data_temp.emplace_back(level,
144 global_dof_indices[i],
145 level_dof_indices[i]);
152 dof_touched[global_idx] =
true;
157 const ::parallel::Triangulation<dim,spacedim> *tria =
159 (&mg_dof.get_triangulation()));
160 AssertThrow(send_data_temp.size()==0 || tria!=
nullptr,
ExcMessage(
"We should only be sending information with a parallel Triangulation!"));
162 #ifdef DEAL_II_WITH_MPI 170 const std::set<types::subdomain_id> &neighbors = tria->level_ghost_owners();
171 std::map<int, std::vector<DoFPair> > send_data;
174 for (
typename std::vector<DoFPair>::iterator dofpair=send_data_temp.begin(); dofpair != send_data_temp.end(); ++dofpair)
176 std::set<types::subdomain_id>::iterator it;
177 for (it = neighbors.begin(); it != neighbors.end(); ++it)
179 if (mg_dof.locally_owned_mg_dofs_per_processor(dofpair->level)[*it].is_element(dofpair->level_dof_index))
181 send_data[*it].push_back(*dofpair);
191 std::vector<MPI_Request> requests;
193 for (std::set<types::subdomain_id>::iterator it = neighbors.begin(); it != neighbors.end(); ++it)
195 requests.push_back(MPI_Request());
196 unsigned int dest = *it;
197 std::vector<DoFPair> &data = send_data[dest];
203 const int ierr = MPI_Isend(data.data(), data.size()*
sizeof(data[0]),
204 MPI_BYTE, dest, 71, tria->get_communicator(),
205 &*requests.rbegin());
210 const int ierr = MPI_Isend(
nullptr, 0, MPI_BYTE, dest, 71,
211 tria->get_communicator(), &*requests.rbegin());
220 std::vector<DoFPair> receive_buffer;
221 for (
unsigned int counter=0; counter<neighbors.size(); ++counter)
225 int ierr = MPI_Probe(MPI_ANY_SOURCE, 71, tria->get_communicator(), &status);
227 ierr = MPI_Get_count(&status, MPI_BYTE, &len);
232 ierr = MPI_Recv(
nullptr, 0, MPI_BYTE, status.MPI_SOURCE, status.MPI_TAG,
233 tria->get_communicator(), &status);
238 int count = len /
sizeof(DoFPair);
240 receive_buffer.resize(count);
242 void *ptr = receive_buffer.data();
243 ierr = MPI_Recv(ptr, len, MPI_BYTE, status.MPI_SOURCE, status.MPI_TAG,
244 tria->get_communicator(), &status);
247 for (
unsigned int i=0; i<receive_buffer.size(); ++i)
249 copy_indices_level_mine[receive_buffer[i].level].emplace_back(
250 receive_buffer[i].global_dof_index, receive_buffer[i].level_dof_index);
256 if (requests.size() > 0)
258 const int ierr = MPI_Waitall(requests.size(), requests.data(), MPI_STATUSES_IGNORE);
266 const int ierr = MPI_Barrier(tria->get_communicator());
275 std::less<std::pair<types::global_dof_index, types::global_dof_index> > compare;
276 for (
unsigned int level=0; level<copy_indices.size(); ++level)
277 std::sort(copy_indices[level].begin(), copy_indices[level].end(), compare);
278 for (
unsigned int level=0; level<copy_indices_level_mine.size(); ++level)
279 std::sort(copy_indices_level_mine[level].begin(), copy_indices_level_mine[level].end(), compare);
280 for (
unsigned int level=0; level<copy_indices_global_mine.size(); ++level)
281 std::sort(copy_indices_global_mine[level].begin(), copy_indices_global_mine[level].end(), compare);
288 template <
typename Number>
290 reinit_ghosted_vector(
const IndexSet &locally_owned,
291 std::vector<types::global_dof_index> &ghosted_level_dofs,
292 const MPI_Comm &communicator,
294 std::vector<std::pair<unsigned int,unsigned int> > ©_indices_global_mine)
296 std::sort(ghosted_level_dofs.begin(), ghosted_level_dofs.end());
298 ghosted_dofs.
add_indices(ghosted_level_dofs.begin(),
299 std::unique(ghosted_level_dofs.begin(),
300 ghosted_level_dofs.end()));
301 ghosted_dofs.compress();
304 if (ghosted_level_vector.
size() == locally_owned.
size())
309 ghosted_dofs.add_indices(part->ghost_indices());
310 for (
unsigned int i=0; i<copy_indices_global_mine.size(); ++i)
311 copy_indices_global_mine[i].second =
313 ghosted_dofs.index_within_set(part->local_to_global(copy_indices_global_mine[i].second));
315 ghosted_level_vector.
reinit(locally_owned, ghosted_dofs, communicator);
321 const std::vector<types::global_dof_index> &mine,
322 const std::vector<types::global_dof_index> &remote,
323 std::vector<unsigned int> &localized_indices)
325 localized_indices.resize(mine.size()+remote.size(),
327 for (
unsigned int i=0; i<mine.size(); ++i)
331 for (
unsigned int i=0; i<remote.size(); ++i)
340 compute_shift_within_children(
const unsigned int child,
341 const unsigned int fe_shift_1d,
342 const unsigned int fe_degree)
346 unsigned int c_tensor_index[dim];
347 unsigned int tmp = child;
348 for (
unsigned int d=0;
d<dim; ++
d)
350 c_tensor_index[
d] = tmp % 2;
353 const unsigned int n_child_dofs_1d = fe_degree + 1 + fe_shift_1d;
354 unsigned int factor = 1;
355 unsigned int shift = fe_shift_1d * c_tensor_index[0];
356 for (
unsigned int d=1;
d<dim; ++
d)
358 factor *= n_child_dofs_1d;
359 shift = shift + factor * fe_shift_1d * c_tensor_index[
d];
367 void add_child_indices(
const unsigned int child,
368 const unsigned int fe_shift_1d,
369 const unsigned int fe_degree,
370 const std::vector<unsigned int> &lexicographic_numbering,
371 const std::vector<types::global_dof_index> &local_dof_indices,
374 const unsigned int n_child_dofs_1d = fe_degree + 1 + fe_shift_1d;
375 const unsigned int shift =
376 compute_shift_within_children<dim>(child, fe_shift_1d, fe_degree);
378 local_dof_indices.size()/Utilities::fixed_power<dim>(fe_degree+1);
380 const unsigned int n_scalar_cell_dofs = Utilities::fixed_power<dim>(n_child_dofs_1d);
382 for (
unsigned int k=0; k<(dim>2 ? (fe_degree+1) : 1); ++k)
383 for (
unsigned int j=0; j<(dim>1 ? (fe_degree+1) : 1); ++j)
384 for (
unsigned int i=0; i<(fe_degree+1); ++i, ++m)
386 const unsigned int index = c*n_scalar_cell_dofs+k*n_child_dofs_1d*
387 n_child_dofs_1d+j*n_child_dofs_1d+i;
389 indices[index] == local_dof_indices[lexicographic_numbering[m]],
391 indices[index] = local_dof_indices[lexicographic_numbering[m]];
395 template <
int dim,
typename Number>
396 void setup_element_info(ElementInfo<Number> &elem_info,
const FiniteElement<1> &fe,
397 const ::DoFHandler<dim> &mg_dof)
400 elem_info.n_components = mg_dof.get_fe().element_multiplicity(0);
402 mg_dof.get_fe().dofs_per_cell);
404 elem_info.fe_degree = fe.
degree;
432 elem_info.n_child_cell_dofs = elem_info.n_components*Utilities::fixed_power<dim>(n_child_dofs_1d);
435 shape_info.
reinit(dummy_quadrature, mg_dof.get_fe(), 0);
439 elem_info.prolongation_matrix_1d.resize(fe.
dofs_per_cell*n_child_dofs_1d);
441 for (
unsigned int c=0; c<GeometryInfo<1>::max_children_per_cell; ++c)
444 elem_info.prolongation_matrix_1d[i*n_child_dofs_1d+j+c*shift] =
455 const unsigned int level,
456 std::vector<types::global_dof_index> &dof_indices)
458 if (mg_constrained_dofs !=
nullptr &&
460 for (
auto &ind : dof_indices)
472 template <
int dim,
typename Number>
473 void setup_transfer(const ::DoFHandler<dim> &mg_dof,
475 ElementInfo<Number> &elem_info,
476 std::vector<std::vector<unsigned int> > &level_dof_indices,
477 std::vector<std::vector<std::pair<unsigned int,unsigned int> > > &parent_child_connect,
478 std::vector<unsigned int> &n_owned_level_cells,
479 std::vector<std::vector<std::vector<unsigned short> > > &dirichlet_indices,
480 std::vector<std::vector<Number> > &weights_on_refined,
481 std::vector<std::vector<std::pair<unsigned int, unsigned int> > > ©_indices_global_mine,
484 level_dof_indices.clear();
485 parent_child_connect.clear();
486 n_owned_level_cells.clear();
487 dirichlet_indices.clear();
488 weights_on_refined.clear();
494 const ::Triangulation<dim> &tria = mg_dof.get_triangulation();
500 std::string fe_name = mg_dof.get_fe().base_element(0).get_name();
502 const std::size_t template_starts = fe_name.find_first_of(
'<');
503 Assert (fe_name[template_starts+1] == (dim==1?
'1':(dim==2?
'2':
'3')),
505 fe_name[template_starts+1] =
'1';
507 const std::unique_ptr<FiniteElement<1> > fe (FETools::get_fe_by_name<1,1>(fe_name));
509 setup_element_info(elem_info,*fe,mg_dof);
513 const unsigned int n_levels = tria.n_global_levels();
514 level_dof_indices.resize(n_levels);
515 parent_child_connect.resize(n_levels-1);
516 n_owned_level_cells.resize(n_levels-1);
517 std::vector<std::vector<unsigned int> > coarse_level_indices(n_levels-1);
518 for (
unsigned int level=0; level<std::min(tria.n_levels(),n_levels-1); ++level)
519 coarse_level_indices[level].resize(tria.n_raw_cells(level),
521 std::vector<types::global_dof_index> local_dof_indices(mg_dof.get_fe().dofs_per_cell);
522 dirichlet_indices.resize(n_levels-1);
527 if (ghosted_level_vector.max_level() != n_levels-1)
528 ghosted_level_vector.resize(0, n_levels-1);
530 for (
unsigned int level=n_levels-1; level > 0; --level)
532 unsigned int counter = 0;
533 std::vector<types::global_dof_index> global_level_dof_indices;
534 std::vector<types::global_dof_index> global_level_dof_indices_remote;
535 std::vector<types::global_dof_index> ghosted_level_dofs;
536 std::vector<types::global_dof_index> global_level_dof_indices_l0;
537 std::vector<types::global_dof_index> ghosted_level_dofs_l0;
540 typename ::DoFHandler<dim>::cell_iterator cell, endc = mg_dof.end(level-1);
541 for (cell = mg_dof.begin(level-1); cell != endc; ++cell)
545 if (!cell->has_children())
548 bool consider_cell =
false;
550 || cell->level_subdomain_id()==tria.locally_owned_subdomain()
552 consider_cell =
true;
557 bool cell_is_remote = !consider_cell;
558 for (
unsigned int c=0; c<GeometryInfo<dim>::max_children_per_cell; ++c)
559 if (cell->child(c)->level_subdomain_id()==tria.locally_owned_subdomain())
561 consider_cell =
true;
576 std::vector<types::global_dof_index> &next_indices =
577 cell_is_remote ? global_level_dof_indices_remote : global_level_dof_indices;
578 const std::size_t start_index = next_indices.size();
579 next_indices.resize(start_index + elem_info.n_child_cell_dofs,
581 for (
unsigned int c=0; c<GeometryInfo<dim>::max_children_per_cell; ++c)
583 if (cell_is_remote && cell->child(c)->level_subdomain_id() !=
584 tria.locally_owned_subdomain())
586 cell->child(c)->get_mg_dof_indices(local_dof_indices);
588 replace(mg_constrained_dofs, level, local_dof_indices);
590 const IndexSet &owned_level_dofs = mg_dof.locally_owned_mg_dofs(level);
591 for (
unsigned int i=0; i<local_dof_indices.size(); ++i)
592 if (!owned_level_dofs.
is_element(local_dof_indices[i]))
593 ghosted_level_dofs.push_back(local_dof_indices[i]);
596 fe->
degree, elem_info.lexicographic_numbering,
598 &next_indices[start_index]);
601 if (cell->child(c)->has_children() &&
603 || cell->child(c)->level_subdomain_id()==tria.locally_owned_subdomain()
606 const unsigned int child_index = coarse_level_indices[level][cell->child(c)->index()];
608 unsigned int parent_index = counter;
616 parent_index = start_index/elem_info.n_child_cell_dofs + tria.n_cells(level);
617 parent_child_connect[level][child_index] =
618 std::make_pair(parent_index, c);
620 static_cast<unsigned short>(-1));
624 if (mg_constrained_dofs !=
nullptr)
625 for (
unsigned int i=0; i<mg_dof.get_fe().dofs_per_cell; ++i)
627 local_dof_indices[elem_info.lexicographic_numbering[i]]))
628 dirichlet_indices[level][child_index].push_back(i);
634 coarse_level_indices[level-1].size());
635 coarse_level_indices[level-1][cell->index()] = counter++;
642 if (level == 1 && !cell_is_remote)
644 cell->get_mg_dof_indices(local_dof_indices);
646 replace(mg_constrained_dofs, level-1, local_dof_indices);
648 const IndexSet &owned_level_dofs_l0 = mg_dof.locally_owned_mg_dofs(0);
649 for (
unsigned int i=0; i<local_dof_indices.size(); ++i)
650 if (!owned_level_dofs_l0.
is_element(local_dof_indices[i]))
651 ghosted_level_dofs_l0.push_back(local_dof_indices[i]);
653 const std::size_t start_index = global_level_dof_indices_l0.size();
654 global_level_dof_indices_l0.resize(start_index+elem_info.n_child_cell_dofs,
657 fe->
degree, elem_info.lexicographic_numbering,
659 &global_level_dof_indices_l0[start_index]);
661 dirichlet_indices[0].emplace_back();
662 if (mg_constrained_dofs !=
nullptr)
663 for (
unsigned int i=0; i<mg_dof.get_fe().dofs_per_cell; ++i)
664 if (mg_constrained_dofs->
is_boundary_index(0, local_dof_indices[elem_info.lexicographic_numbering[i]]))
665 dirichlet_indices[0].back().push_back(i);
672 AssertDimension(counter*elem_info.n_child_cell_dofs, global_level_dof_indices.size());
673 n_owned_level_cells[level-1] = counter;
674 dirichlet_indices[level-1].resize(counter);
675 parent_child_connect[level-1].
682 if (level < n_levels-1)
683 for (std::vector<std::pair<unsigned int,unsigned int> >::iterator
684 i=parent_child_connect[level].begin(); i!=parent_child_connect[level].end(); ++i)
685 if (i->first >= tria.n_cells(level))
687 i->first -= tria.n_cells(level);
694 const MPI_Comm communicator =
697 reinit_ghosted_vector (mg_dof.locally_owned_mg_dofs(level),
698 ghosted_level_dofs, communicator,
699 ghosted_level_vector[level],
700 copy_indices_global_mine[level]);
702 copy_indices_to_mpi_local_numbers(*ghosted_level_vector[level].get_partitioner(),
703 global_level_dof_indices,
704 global_level_dof_indices_remote,
705 level_dof_indices[level]);
709 for (
unsigned int i = 0; i<parent_child_connect[0].size(); ++i)
710 parent_child_connect[0][i] = std::make_pair(i, 0U);
712 reinit_ghosted_vector (mg_dof.locally_owned_mg_dofs(0),
713 ghosted_level_dofs_l0, communicator,
714 ghosted_level_vector[0],
715 copy_indices_global_mine[0]);
717 copy_indices_to_mpi_local_numbers(*ghosted_level_vector[0].get_partitioner(),
718 global_level_dof_indices_l0,
719 std::vector<types::global_dof_index>(),
720 level_dof_indices[0]);
730 weights_on_refined.resize(n_levels-1);
731 for (
unsigned int level = 1; level<n_levels; ++level)
733 ghosted_level_vector[level] = 0;
734 for (
unsigned int c=0; c<n_owned_level_cells[level-1]; ++c)
735 for (
unsigned int j=0; j<elem_info.n_child_cell_dofs; ++j)
736 ghosted_level_vector[level].local_element(level_dof_indices[level][elem_info.n_child_cell_dofs*c+j]) += Number(1.);
738 ghosted_level_vector[level].update_ghost_values();
740 std::vector<unsigned int> degree_to_3 (n_child_dofs_1d);
742 for (
unsigned int i=1; i<n_child_dofs_1d-1; ++i)
744 degree_to_3.back() = 2;
748 weights_on_refined[level-1].resize(n_owned_level_cells[level-1]*Utilities::fixed_power<dim>(3));
749 for (
unsigned int c=0; c<n_owned_level_cells[level-1]; ++c)
750 for (
unsigned int k=0, m=0; k<(dim>2 ? n_child_dofs_1d : 1); ++k)
751 for (
unsigned int j=0; j<(dim>1 ? n_child_dofs_1d : 1); ++j)
753 unsigned int shift = 9*degree_to_3[k] + 3*degree_to_3[j];
754 for (
unsigned int i=0; i<n_child_dofs_1d; ++i, ++m)
755 weights_on_refined[level-1][c*Utilities::fixed_power<dim>(3)+shift+degree_to_3[i]] = Number(1.)/
756 ghosted_level_vector[level].local_element(level_dof_indices[level][elem_info.n_child_cell_dofs*c+m]);
767 #include "mg_transfer_internal.inst" 769 DEAL_II_NAMESPACE_CLOSE
void reinit(const size_type size, const bool omit_zeroing_entries=false)
static const unsigned int invalid_unsigned_int
#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)
bool is_boundary_index(const unsigned int level, const types::global_dof_index index) const
#define AssertIndexRange(index, range)
void add_indices(const ForwardIterator &begin, const ForwardIterator &end)
const unsigned int degree
#define AssertThrow(cond, exc)
virtual const FullMatrix< double > & get_prolongation_matrix(const unsigned int child, const RefinementCase< dim > &refinement_case=RefinementCase< dim >::isotropic_refinement) const
unsigned int global_to_local(const types::global_dof_index global_index) const
const unsigned int dofs_per_line
static ::ExceptionBase & ExcMessage(std::string arg1)
unsigned int global_dof_index
virtual size_type size() const override
#define Assert(cond, exc)
size_type index_within_set(const size_type global_index) const
size_type n_constraints() const
virtual MPI_Comm get_communicator() const
const ConstraintMatrix & get_level_constraint_matrix(const unsigned int level) const
const std::shared_ptr< const Utilities::MPI::Partitioner > & get_partitioner() const
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
const unsigned int dofs_per_cell
const std::vector< Point< dim > > & get_unit_support_points() const
const types::subdomain_id artificial_subdomain_id
bool at_refinement_edge(const unsigned int level, const types::global_dof_index index) const
#define AssertThrowMPI(error_code)
bool is_identity_constrained(const size_type index) const
static ::ExceptionBase & ExcNotImplemented()
const unsigned int dofs_per_vertex
bool is_element(const size_type index) const
const std::vector< std::pair< size_type, double > > * get_constraint_entries(const size_type line) const
std::vector< unsigned int > lexicographic_numbering
const types::global_dof_index invalid_dof_index
size_type n_elements() const
static ::ExceptionBase & ExcInternalError()