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),
58 template <
int dim,
int spacedim>
59 void fill_copy_indices(const ::DoFHandler<dim,spacedim> &mg_dof,
61 std::vector<std::vector<std::pair<types::global_dof_index, types::global_dof_index> > > ©_indices,
62 std::vector<std::vector<std::pair<types::global_dof_index, types::global_dof_index> > > ©_indices_global_mine,
63 std::vector<std::vector<std::pair<types::global_dof_index, types::global_dof_index> > > ©_indices_level_mine)
76 std::vector<DoFPair> send_data_temp;
78 const unsigned int n_levels = mg_dof.get_triangulation().n_global_levels();
79 copy_indices.resize(n_levels);
80 copy_indices_global_mine.resize(n_levels);
81 copy_indices_level_mine.resize(n_levels);
85 const unsigned int dofs_per_cell = mg_dof.get_fe().dofs_per_cell;
86 std::vector<types::global_dof_index> global_dof_indices (dofs_per_cell);
87 std::vector<types::global_dof_index> level_dof_indices (dofs_per_cell);
89 for (
unsigned int level=0; level<n_levels; ++level)
91 std::vector<bool> dof_touched(globally_relevant.
n_elements(),
false);
92 copy_indices[level].clear();
93 copy_indices_level_mine[level].clear();
94 copy_indices_global_mine[level].clear();
96 typename ::DoFHandler<dim,spacedim>::active_cell_iterator
97 level_cell = mg_dof.begin_active(level);
98 const typename ::DoFHandler<dim,spacedim>::active_cell_iterator
99 level_end = mg_dof.end_active(level);
101 for (; level_cell!=level_end; ++level_cell)
111 level_cell->get_dof_indices (global_dof_indices);
112 level_cell->get_mg_dof_indices (level_dof_indices);
114 for (
unsigned int i=0; i<dofs_per_cell; ++i)
117 if (mg_constrained_dofs != 0
122 if (dof_touched[global_idx])
124 bool global_mine = mg_dof.locally_owned_dofs().is_element(global_dof_indices[i]);
125 bool level_mine = mg_dof.locally_owned_mg_dofs(level).is_element(level_dof_indices[i]);
128 if (global_mine && level_mine)
130 copy_indices[level].push_back(
131 std::make_pair (global_dof_indices[i], level_dof_indices[i]));
133 else if (global_mine)
135 copy_indices_global_mine[level].push_back(
136 std::make_pair (global_dof_indices[i], level_dof_indices[i]));
139 send_data_temp.push_back(DoFPair(level, global_dof_indices[i], level_dof_indices[i]));
146 dof_touched[global_idx] =
true;
151 const ::parallel::distributed::Triangulation<dim,spacedim> *tria =
153 (&mg_dof.get_triangulation()));
154 AssertThrow(send_data_temp.size()==0 || tria!=NULL,
ExcMessage(
"parallel Multigrid only works with a distributed Triangulation!"));
156 #ifdef DEAL_II_WITH_MPI 164 std::set<types::subdomain_id> neighbors = tria->level_ghost_owners();
165 std::map<int, std::vector<DoFPair> > send_data;
168 for (
typename std::vector<DoFPair>::iterator dofpair=send_data_temp.begin(); dofpair != send_data_temp.end(); ++dofpair)
170 std::set<types::subdomain_id>::iterator it;
171 for (it = neighbors.begin(); it != neighbors.end(); ++it)
173 if (mg_dof.locally_owned_mg_dofs_per_processor(dofpair->level)[*it].is_element(dofpair->level_dof_index))
175 send_data[*it].push_back(*dofpair);
185 std::vector<MPI_Request> requests;
187 for (std::set<types::subdomain_id>::iterator it = neighbors.begin(); it != neighbors.end(); ++it)
189 requests.push_back(MPI_Request());
190 unsigned int dest = *it;
191 std::vector<DoFPair> &data = send_data[dest];
197 const int ierr = MPI_Isend(&data[0], data.size()*
sizeof(data[0]),
198 MPI_BYTE, dest, 71, tria->get_communicator(),
199 &*requests.rbegin());
204 const int ierr = MPI_Isend(NULL, 0, MPI_BYTE, dest, 71,
205 tria->get_communicator(), &*requests.rbegin());
214 std::vector<DoFPair> receive_buffer;
215 for (
unsigned int counter=0; counter<neighbors.size(); ++counter)
219 int ierr = MPI_Probe(MPI_ANY_SOURCE, 71, tria->get_communicator(), &status);
221 ierr = MPI_Get_count(&status, MPI_BYTE, &len);
226 ierr = MPI_Recv(NULL, 0, MPI_BYTE, status.MPI_SOURCE, status.MPI_TAG,
227 tria->get_communicator(), &status);
232 int count = len /
sizeof(DoFPair);
234 receive_buffer.resize(count);
236 void *ptr = &receive_buffer[0];
237 ierr = MPI_Recv(ptr, len, MPI_BYTE, status.MPI_SOURCE, status.MPI_TAG,
238 tria->get_communicator(), &status);
241 for (
unsigned int i=0; i<receive_buffer.size(); ++i)
243 copy_indices_level_mine[receive_buffer[i].level].push_back(
244 std::make_pair (receive_buffer[i].global_dof_index, receive_buffer[i].level_dof_index)
251 if (requests.size() > 0)
253 const int ierr = MPI_Waitall(requests.size(), &requests[0], MPI_STATUSES_IGNORE);
261 const int ierr = MPI_Barrier(tria->get_communicator());
270 std::less<std::pair<types::global_dof_index, types::global_dof_index> > compare;
271 for (
unsigned int level=0; level<copy_indices.size(); ++level)
272 std::sort(copy_indices[level].begin(), copy_indices[level].end(), compare);
273 for (
unsigned int level=0; level<copy_indices_level_mine.size(); ++level)
274 std::sort(copy_indices_level_mine[level].begin(), copy_indices_level_mine[level].end(), compare);
275 for (
unsigned int level=0; level<copy_indices_global_mine.size(); ++level)
276 std::sort(copy_indices_global_mine[level].begin(), copy_indices_global_mine[level].end(), compare);
283 template <
typename Number>
285 reinit_ghosted_vector(
const IndexSet &locally_owned,
286 std::vector<types::global_dof_index> &ghosted_level_dofs,
287 const MPI_Comm &communicator,
289 std::vector<std::pair<unsigned int,unsigned int> > ©_indices_global_mine)
291 std::sort(ghosted_level_dofs.begin(), ghosted_level_dofs.end());
293 ghosted_dofs.
add_indices(ghosted_level_dofs.begin(),
294 std::unique(ghosted_level_dofs.begin(),
295 ghosted_level_dofs.end()));
296 ghosted_dofs.compress();
299 if (ghosted_level_vector.
size() == locally_owned.
size())
303 const std_cxx11::shared_ptr<const Utilities::MPI::Partitioner> part
305 ghosted_dofs.add_indices(part->ghost_indices());
306 for (
unsigned int i=0; i<copy_indices_global_mine.size(); ++i)
307 copy_indices_global_mine[i].second =
309 ghosted_dofs.index_within_set(part->local_to_global(copy_indices_global_mine[i].second));
311 ghosted_level_vector.
reinit(locally_owned, ghosted_dofs, communicator);
317 const std::vector<types::global_dof_index> &mine,
318 const std::vector<types::global_dof_index> &remote,
319 std::vector<unsigned int> &localized_indices)
321 localized_indices.resize(mine.size()+remote.size(),
323 for (
unsigned int i=0; i<mine.size(); ++i)
327 for (
unsigned int i=0; i<remote.size(); ++i)
336 compute_shift_within_children(
const unsigned int child,
337 const unsigned int fe_shift_1d,
338 const unsigned int fe_degree)
342 unsigned int c_tensor_index[dim];
343 unsigned int tmp = child;
344 for (
unsigned int d=0;
d<dim; ++
d)
346 c_tensor_index[
d] = tmp % 2;
349 const unsigned int n_child_dofs_1d = fe_degree + 1 + fe_shift_1d;
350 unsigned int factor = 1;
351 unsigned int shift = fe_shift_1d * c_tensor_index[0];
352 for (
unsigned int d=1;
d<dim; ++
d)
354 factor *= n_child_dofs_1d;
355 shift = shift + factor * fe_shift_1d * c_tensor_index[
d];
363 void add_child_indices(
const unsigned int child,
364 const unsigned int fe_shift_1d,
365 const unsigned int fe_degree,
366 const std::vector<unsigned int> &lexicographic_numbering,
367 const std::vector<types::global_dof_index> &local_dof_indices,
370 const unsigned int n_child_dofs_1d = fe_degree + 1 + fe_shift_1d;
371 const unsigned int shift =
372 compute_shift_within_children<dim>(child, fe_shift_1d, fe_degree);
373 const unsigned int n_components =
374 local_dof_indices.size()/Utilities::fixed_power<dim>(fe_degree+1);
376 const unsigned int n_scalar_cell_dofs = Utilities::fixed_power<dim>(n_child_dofs_1d);
377 for (
unsigned int c=0, m=0; c<n_components; ++c)
378 for (
unsigned int k=0; k<(dim>2 ? (fe_degree+1) : 1); ++k)
379 for (
unsigned int j=0; j<(dim>1 ? (fe_degree+1) : 1); ++j)
380 for (
unsigned int i=0; i<(fe_degree+1); ++i, ++m)
382 const unsigned int index = c*n_scalar_cell_dofs+k*n_child_dofs_1d*
383 n_child_dofs_1d+j*n_child_dofs_1d+i;
385 indices[index] == local_dof_indices[lexicographic_numbering[m]],
387 indices[index] = local_dof_indices[lexicographic_numbering[m]];
391 template <
int dim,
typename Number>
392 void setup_element_info(ElementInfo<Number> &elem_info,
const FiniteElement<1> &fe,
393 const ::DoFHandler<dim> &mg_dof)
396 elem_info.n_components = mg_dof.get_fe().element_multiplicity(0);
398 mg_dof.get_fe().dofs_per_cell);
400 elem_info.fe_degree = fe.
degree;
428 elem_info.n_child_cell_dofs = elem_info.n_components*Utilities::fixed_power<dim>(n_child_dofs_1d);
431 shape_info.
reinit(dummy_quadrature, mg_dof.get_fe(), 0);
435 elem_info.prolongation_matrix_1d.resize(fe.
dofs_per_cell*n_child_dofs_1d);
437 for (
unsigned int c=0; c<GeometryInfo<1>::max_children_per_cell; ++c)
440 elem_info.prolongation_matrix_1d[i*n_child_dofs_1d+j+c*shift] =
447 template <
int dim,
typename Number>
448 void setup_transfer(const ::DoFHandler<dim> &mg_dof,
450 ElementInfo<Number> &elem_info,
451 std::vector<std::vector<unsigned int> > &level_dof_indices,
452 std::vector<std::vector<std::pair<unsigned int,unsigned int> > > &parent_child_connect,
453 std::vector<unsigned int> &n_owned_level_cells,
454 std::vector<std::vector<std::vector<unsigned short> > > &dirichlet_indices,
455 std::vector<std::vector<Number> > &weights_on_refined,
456 std::vector<std::vector<std::pair<unsigned int, unsigned int> > > ©_indices_global_mine,
459 level_dof_indices.clear();
460 parent_child_connect.clear();
461 n_owned_level_cells.clear();
462 dirichlet_indices.clear();
463 weights_on_refined.clear();
469 const ::Triangulation<dim> &tria = mg_dof.get_triangulation();
475 std::string fe_name = mg_dof.get_fe().base_element(0).get_name();
477 const std::size_t template_starts = fe_name.find_first_of(
'<');
478 Assert (fe_name[template_starts+1] == (dim==1?
'1':(dim==2?
'2':
'3')),
480 fe_name[template_starts+1] =
'1';
482 std_cxx11::shared_ptr<FiniteElement<1> > fe_1d
483 (FETools::get_fe_by_name<1,1>(fe_name));
486 setup_element_info(elem_info,fe,mg_dof);
490 const unsigned int n_levels = tria.n_global_levels();
491 level_dof_indices.resize(n_levels);
492 parent_child_connect.resize(n_levels-1);
493 n_owned_level_cells.resize(n_levels-1);
494 std::vector<std::vector<unsigned int> > coarse_level_indices(n_levels-1);
495 for (
unsigned int level=0; level<std::min(tria.n_levels(),n_levels-1); ++level)
496 coarse_level_indices[level].resize(tria.n_raw_cells(level),
498 std::vector<types::global_dof_index> local_dof_indices(mg_dof.get_fe().dofs_per_cell);
499 dirichlet_indices.resize(n_levels-1);
504 if (ghosted_level_vector.max_level() != n_levels-1)
505 ghosted_level_vector.resize(0, n_levels-1);
507 for (
unsigned int level=n_levels-1; level > 0; --level)
509 unsigned int counter = 0;
510 std::vector<types::global_dof_index> global_level_dof_indices;
511 std::vector<types::global_dof_index> global_level_dof_indices_remote;
512 std::vector<types::global_dof_index> ghosted_level_dofs;
513 std::vector<types::global_dof_index> global_level_dof_indices_l0;
514 std::vector<types::global_dof_index> ghosted_level_dofs_l0;
517 typename ::DoFHandler<dim>::cell_iterator cell, endc = mg_dof.end(level-1);
518 for (cell = mg_dof.begin(level-1); cell != endc; ++cell)
522 if (!cell->has_children())
525 bool consider_cell =
false;
527 || cell->level_subdomain_id()==tria.locally_owned_subdomain()
529 consider_cell =
true;
534 bool cell_is_remote = !consider_cell;
535 for (
unsigned int c=0; c<GeometryInfo<dim>::max_children_per_cell; ++c)
536 if (cell->child(c)->level_subdomain_id()==tria.locally_owned_subdomain())
538 consider_cell =
true;
553 std::vector<types::global_dof_index> &next_indices =
554 cell_is_remote ? global_level_dof_indices_remote : global_level_dof_indices;
555 const std::size_t start_index = next_indices.size();
556 next_indices.resize(start_index + elem_info.n_child_cell_dofs,
558 for (
unsigned int c=0; c<GeometryInfo<dim>::max_children_per_cell; ++c)
560 if (cell_is_remote && cell->child(c)->level_subdomain_id() !=
561 tria.locally_owned_subdomain())
563 cell->child(c)->get_mg_dof_indices(local_dof_indices);
565 const IndexSet &owned_level_dofs = mg_dof.locally_owned_mg_dofs(level);
566 for (
unsigned int i=0; i<local_dof_indices.size(); ++i)
567 if (!owned_level_dofs.
is_element(local_dof_indices[i]))
568 ghosted_level_dofs.push_back(local_dof_indices[i]);
571 fe.
degree, elem_info.lexicographic_numbering,
573 &next_indices[start_index]);
576 if (cell->child(c)->has_children() &&
578 || cell->child(c)->level_subdomain_id()==tria.locally_owned_subdomain()
581 const unsigned int child_index = coarse_level_indices[level][cell->child(c)->index()];
583 unsigned int parent_index = counter;
591 parent_index = start_index/elem_info.n_child_cell_dofs + tria.n_cells(level);
592 parent_child_connect[level][child_index] =
593 std::make_pair(parent_index, c);
595 static_cast<unsigned short>(-1));
599 if (mg_constrained_dofs != 0)
600 for (
unsigned int i=0; i<mg_dof.get_fe().dofs_per_cell; ++i)
602 local_dof_indices[elem_info.lexicographic_numbering[i]]))
603 dirichlet_indices[level][child_index].push_back(i);
609 coarse_level_indices[level-1].size());
610 coarse_level_indices[level-1][cell->index()] = counter++;
617 if (level == 1 && !cell_is_remote)
619 cell->get_mg_dof_indices(local_dof_indices);
621 const IndexSet &owned_level_dofs_l0 = mg_dof.locally_owned_mg_dofs(0);
622 for (
unsigned int i=0; i<local_dof_indices.size(); ++i)
623 if (!owned_level_dofs_l0.
is_element(local_dof_indices[i]))
624 ghosted_level_dofs_l0.push_back(local_dof_indices[i]);
626 const std::size_t start_index = global_level_dof_indices_l0.size();
627 global_level_dof_indices_l0.resize(start_index+elem_info.n_child_cell_dofs,
630 fe.
degree, elem_info.lexicographic_numbering,
632 &global_level_dof_indices_l0[start_index]);
634 dirichlet_indices[0].push_back(std::vector<unsigned short>());
635 if (mg_constrained_dofs != 0)
636 for (
unsigned int i=0; i<mg_dof.get_fe().dofs_per_cell; ++i)
637 if (mg_constrained_dofs->
is_boundary_index(0, local_dof_indices[elem_info.lexicographic_numbering[i]]))
638 dirichlet_indices[0].back().push_back(i);
645 AssertDimension(counter*elem_info.n_child_cell_dofs, global_level_dof_indices.size());
646 n_owned_level_cells[level-1] = counter;
647 dirichlet_indices[level-1].resize(counter);
648 parent_child_connect[level-1].
655 if (level < n_levels-1)
656 for (std::vector<std::pair<unsigned int,unsigned int> >::iterator
657 i=parent_child_connect[level].begin(); i!=parent_child_connect[level].end(); ++i)
658 if (i->first >= tria.n_cells(level))
660 i->first -= tria.n_cells(level);
667 const MPI_Comm communicator =
670 reinit_ghosted_vector (mg_dof.locally_owned_mg_dofs(level),
671 ghosted_level_dofs, communicator,
672 ghosted_level_vector[level],
673 copy_indices_global_mine[level]);
675 copy_indices_to_mpi_local_numbers(*ghosted_level_vector[level].get_partitioner(),
676 global_level_dof_indices,
677 global_level_dof_indices_remote,
678 level_dof_indices[level]);
682 for (
unsigned int i = 0; i<parent_child_connect[0].size(); ++i)
683 parent_child_connect[0][i] = std::make_pair(i, 0U);
685 reinit_ghosted_vector (mg_dof.locally_owned_mg_dofs(0),
686 ghosted_level_dofs_l0, communicator,
687 ghosted_level_vector[0],
688 copy_indices_global_mine[0]);
690 copy_indices_to_mpi_local_numbers(*ghosted_level_vector[0].get_partitioner(),
691 global_level_dof_indices_l0,
692 std::vector<types::global_dof_index>(),
693 level_dof_indices[0]);
703 weights_on_refined.resize(n_levels-1);
704 for (
unsigned int level = 1; level<n_levels; ++level)
706 ghosted_level_vector[level] = 0;
707 for (
unsigned int c=0; c<n_owned_level_cells[level-1]; ++c)
708 for (
unsigned int j=0; j<elem_info.n_child_cell_dofs; ++j)
709 ghosted_level_vector[level].local_element(level_dof_indices[level][elem_info.n_child_cell_dofs*c+j]) += Number(1.);
711 ghosted_level_vector[level].update_ghost_values();
713 std::vector<unsigned int> degree_to_3 (n_child_dofs_1d);
715 for (
unsigned int i=1; i<n_child_dofs_1d-1; ++i)
717 degree_to_3.back() = 2;
721 weights_on_refined[level-1].resize(n_owned_level_cells[level-1]*Utilities::fixed_power<dim>(3));
722 for (
unsigned int c=0; c<n_owned_level_cells[level-1]; ++c)
723 for (
unsigned int k=0, m=0; k<(dim>2 ? n_child_dofs_1d : 1); ++k)
724 for (
unsigned int j=0; j<(dim>1 ? n_child_dofs_1d : 1); ++j)
726 unsigned int shift = 9*degree_to_3[k] + 3*degree_to_3[j];
727 for (
unsigned int i=0; i<n_child_dofs_1d; ++i, ++m)
728 weights_on_refined[level-1][c*Utilities::fixed_power<dim>(3)+shift+degree_to_3[i]] = Number(1.)/
729 ghosted_level_vector[level].local_element(level_dof_indices[level][elem_info.n_child_cell_dofs*c+m]);
740 #include "mg_transfer_internal.inst" 742 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)
virtual size_type size() const
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
#define Assert(cond, exc)
size_type index_within_set(const size_type global_index) const
virtual MPI_Comm get_communicator() 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)
static ::ExceptionBase & ExcNotImplemented()
const unsigned int dofs_per_vertex
bool is_element(const size_type index) const
std::vector< unsigned int > lexicographic_numbering
const types::global_dof_index invalid_dof_index
const std_cxx11::shared_ptr< const Utilities::MPI::Partitioner > & get_partitioner() const
size_type n_elements() const
static ::ExceptionBase & ExcInternalError()