36 #include <deal.II/multigrid/mg_transfer_matrix_free.templates.h>
43 template <
int dim,
typename Number>
46 , element_is_continuous(false)
48 , n_child_cell_dofs(0)
53 template <
int dim,
typename Number>
57 , element_is_continuous(false)
59 , n_child_cell_dofs(0)
66 template <
int dim,
typename Number>
71 this->mg_constrained_dofs = &mg_c;
76 template <
int dim,
typename Number>
83 element_is_continuous =
false;
85 n_child_cell_dofs = 0;
86 level_dof_indices.clear();
87 parent_child_connect.clear();
88 dirichlet_indices.clear();
89 n_owned_level_cells.clear();
90 prolongation_matrix_1d.clear();
91 evaluation_data.clear();
92 weights_on_refined.clear();
97 template <
int dim,
typename Number>
101 const std::vector<std::shared_ptr<const Utilities::MPI::Partitioner>>
102 &external_partitioners)
106 "The underlying DoFHandler object has not had its "
107 "distribute_mg_dofs() function called, but this is a prerequisite "
108 "for multigrid transfers. You will need to call this function, "
109 "probably close to where you already call distribute_dofs()."));
110 if (external_partitioners.size() > 0)
113 this->initialize_dof_vector ==
nullptr,
115 "A initialize_dof_vector function has already been registered in the constructor!"));
117 this->initialize_dof_vector =
118 [external_partitioners](
119 const unsigned int level,
121 vec.reinit(external_partitioners[
level]);
125 this->fill_and_communicate_copy_indices(dof_handler);
127 vector_partitioners.resize(0,
130 for (
unsigned int level = 0;
level <= this->ghosted_level_vector.max_level();
132 vector_partitioners[
level] =
133 this->ghosted_level_vector[
level].get_partitioner();
135 std::vector<std::vector<Number>> weights_unvectorized;
139 internal::MGTransfer::setup_transfer<dim, Number>(
141 this->mg_constrained_dofs,
142 external_partitioners,
145 parent_child_connect,
148 weights_unvectorized,
149 this->copy_indices_global_mine,
150 vector_partitioners);
154 this->ghosted_level_vector.resize(0, vector_partitioners.max_level());
155 for (
unsigned int level = 0;
level <= vector_partitioners.max_level();
157 if (external_partitioners.size() == vector_partitioners.max_level() + 1 &&
158 external_partitioners[
level].get() == vector_partitioners[
level].get())
159 this->ghosted_level_vector[
level].reinit(0);
161 this->ghosted_level_vector[
level].reinit(vector_partitioners[
level]);
176 const unsigned int n_levels =
179 const unsigned int n_weights_per_cell = Utilities::fixed_power<dim>(3);
180 weights_on_refined.resize(n_levels - 1);
183 weights_on_refined[
level - 1].resize(
184 ((n_owned_level_cells[
level - 1] + vec_size - 1) / vec_size) *
187 for (
unsigned int c = 0; c < n_owned_level_cells[
level - 1]; ++c)
189 const unsigned int comp = c / vec_size;
190 const unsigned int v = c % vec_size;
191 for (
unsigned int i = 0; i < n_weights_per_cell; ++i)
193 weights_on_refined[
level - 1][comp * n_weights_per_cell + i][v] =
194 weights_unvectorized[
level - 1][c * n_weights_per_cell + i];
199 evaluation_data.resize(n_child_cell_dofs);
204 template <
int dim,
typename Number>
208 const std::function<
void(
const unsigned int,
210 &initialize_dof_vector)
212 if (initialize_dof_vector)
214 const unsigned int n_levels =
217 std::vector<std::shared_ptr<const Utilities::MPI::Partitioner>>
218 external_partitioners(n_levels);
223 initialize_dof_vector(
level, vector);
227 build(dof_handler, external_partitioners);
237 template <
int dim,
typename Number>
240 const unsigned int to_level,
245 prolongate_and_add(to_level, dst, src);
250 template <
int dim,
typename Number>
253 const unsigned int to_level,
257 Assert((to_level >= 1) && (to_level <= level_dof_indices.size()),
261 this->vector_partitioners[to_level - 1].get();
262 if (src_inplace ==
false)
264 if (this->ghosted_level_vector[to_level - 1].get_partitioner().get() !=
265 this->vector_partitioners[to_level - 1].get())
266 this->ghosted_level_vector[to_level - 1].reinit(
267 this->vector_partitioners[to_level - 1]);
268 this->ghosted_level_vector[to_level - 1].copy_locally_owned_data_from(
272 const bool dst_inplace =
273 dst.
get_partitioner().get() == this->vector_partitioners[to_level].get();
274 if (dst_inplace ==
false)
276 if (this->ghosted_level_vector[to_level].get_partitioner().get() !=
277 this->vector_partitioners[to_level].get())
278 this->ghosted_level_vector[to_level].reinit(
279 this->vector_partitioners[to_level]);
280 AssertDimension(this->ghosted_level_vector[to_level].locally_owned_size(),
282 this->ghosted_level_vector[to_level] = 0.;
286 src_inplace ? src : this->ghosted_level_vector[to_level - 1];
288 dst_inplace ? dst : this->ghosted_level_vector[to_level];
295 do_prolongate_add<0>(to_level, dst_vec, src_vec);
296 else if (fe_degree == 1)
297 do_prolongate_add<1>(to_level, dst_vec, src_vec);
298 else if (fe_degree == 2)
299 do_prolongate_add<2>(to_level, dst_vec, src_vec);
300 else if (fe_degree == 3)
301 do_prolongate_add<3>(to_level, dst_vec, src_vec);
302 else if (fe_degree == 4)
303 do_prolongate_add<4>(to_level, dst_vec, src_vec);
304 else if (fe_degree == 5)
305 do_prolongate_add<5>(to_level, dst_vec, src_vec);
306 else if (fe_degree == 6)
307 do_prolongate_add<6>(to_level, dst_vec, src_vec);
308 else if (fe_degree == 7)
309 do_prolongate_add<7>(to_level, dst_vec, src_vec);
310 else if (fe_degree == 8)
311 do_prolongate_add<8>(to_level, dst_vec, src_vec);
312 else if (fe_degree == 9)
313 do_prolongate_add<9>(to_level, dst_vec, src_vec);
314 else if (fe_degree == 10)
315 do_prolongate_add<10>(to_level, dst_vec, src_vec);
317 do_prolongate_add<-1>(to_level, dst_vec, src_vec);
320 if (dst_inplace ==
false)
323 if (src_inplace ==
true)
329 template <
int dim,
typename Number>
332 const unsigned int from_level,
336 Assert((from_level >= 1) && (from_level <= level_dof_indices.size()),
339 const bool src_inplace =
340 src.
get_partitioner().get() == this->vector_partitioners[from_level].get();
341 if (src_inplace ==
false)
343 if (this->ghosted_level_vector[from_level].get_partitioner().get() !=
344 this->vector_partitioners[from_level].get())
345 this->ghosted_level_vector[from_level].reinit(
346 this->vector_partitioners[from_level]);
347 this->ghosted_level_vector[from_level].copy_locally_owned_data_from(src);
351 this->vector_partitioners[from_level - 1].get();
352 if (dst_inplace ==
false)
354 if (this->ghosted_level_vector[from_level - 1].get_partitioner().get() !=
355 this->vector_partitioners[from_level - 1].get())
356 this->ghosted_level_vector[from_level - 1].reinit(
357 this->vector_partitioners[from_level - 1]);
359 this->ghosted_level_vector[from_level - 1].locally_owned_size(),
361 this->ghosted_level_vector[from_level - 1] = 0.;
365 src_inplace ? src : this->ghosted_level_vector[from_level];
367 dst_inplace ? dst : this->ghosted_level_vector[from_level - 1];
372 do_restrict_add<0>(from_level, dst_vec, src_vec);
373 else if (fe_degree == 1)
374 do_restrict_add<1>(from_level, dst_vec, src_vec);
375 else if (fe_degree == 2)
376 do_restrict_add<2>(from_level, dst_vec, src_vec);
377 else if (fe_degree == 3)
378 do_restrict_add<3>(from_level, dst_vec, src_vec);
379 else if (fe_degree == 4)
380 do_restrict_add<4>(from_level, dst_vec, src_vec);
381 else if (fe_degree == 5)
382 do_restrict_add<5>(from_level, dst_vec, src_vec);
383 else if (fe_degree == 6)
384 do_restrict_add<6>(from_level, dst_vec, src_vec);
385 else if (fe_degree == 7)
386 do_restrict_add<7>(from_level, dst_vec, src_vec);
387 else if (fe_degree == 8)
388 do_restrict_add<8>(from_level, dst_vec, src_vec);
389 else if (fe_degree == 9)
390 do_restrict_add<9>(from_level, dst_vec, src_vec);
391 else if (fe_degree == 10)
392 do_restrict_add<10>(from_level, dst_vec, src_vec);
395 do_restrict_add<-1>(from_level, dst_vec, src_vec);
398 if (dst_inplace ==
false)
401 if (src_inplace ==
true)
407 template <
int dim,
typename Number>
408 template <
int degree>
411 const unsigned int to_level,
416 const unsigned int degree_size = (degree > -1 ? degree : fe_degree) + 1;
417 const unsigned int n_child_dofs_1d = 2 * degree_size - element_is_continuous;
418 const unsigned int n_scalar_cell_dofs =
419 Utilities::fixed_power<dim>(n_child_dofs_1d);
427 if (this->mg_constrained_dofs !=
nullptr &&
428 this->mg_constrained_dofs->get_user_constraint_matrix(to_level - 1)
435 this->mg_constrained_dofs->get_user_constraint_matrix(to_level - 1)
436 .distribute(copy_src);
446 for (
unsigned int cell = 0; cell < n_owned_level_cells[to_level - 1];
449 const unsigned int n_lanes =
450 (cell + vec_size > n_owned_level_cells[to_level - 1]) ?
451 (n_owned_level_cells[to_level - 1] - cell) :
455 for (
unsigned int v = 0; v < n_lanes; ++v)
457 const unsigned int shift =
458 internal::MGTransfer::compute_shift_within_children<dim>(
459 parent_child_connect[to_level - 1][cell + v].
second,
460 fe_degree + 1 - element_is_continuous,
462 const unsigned int *indices =
463 &level_dof_indices[to_level - 1]
464 [parent_child_connect[to_level - 1][cell + v]
468 for (
unsigned int c = 0, m = 0; c < n_components; ++c)
470 for (
unsigned int k = 0; k < (dim > 2 ? degree_size : 1); ++k)
471 for (
unsigned int j = 0; j < (dim > 1 ? degree_size : 1); ++j)
472 for (
unsigned int i = 0; i < degree_size; ++i, ++m)
474 indices[c * n_scalar_cell_dofs +
475 k * n_child_dofs_1d * n_child_dofs_1d +
476 j * n_child_dofs_1d + i]);
479 for (std::vector<unsigned short>::const_iterator i =
480 dirichlet_indices[to_level - 1][cell + v].
begin();
481 i != dirichlet_indices[to_level - 1][cell + v].end();
483 evaluation_data[*i][v] = 0.;
488 degree_size * n_child_dofs_1d);
490 if (element_is_continuous)
494 for (
int c = n_components - 1; c >= 0; --c)
503 prolongation_matrix_1d,
504 evaluation_data.begin() +
507 evaluation_data.begin() +
508 c * n_scalar_cell_dofs,
512 degree != -1 ? 2 * degree + 1 :
515 &weights_on_refined[to_level - 1][(cell / vec_size) * three_to_dim],
518 evaluation_data.begin());
522 for (
int c = n_components - 1; c >= 0; --c)
531 prolongation_matrix_1d,
532 evaluation_data.begin() +
535 evaluation_data.begin() +
536 c * n_scalar_cell_dofs,
542 const unsigned int *indices =
543 &level_dof_indices[to_level][cell * n_child_cell_dofs];
544 for (
unsigned int v = 0; v < n_lanes; ++v)
546 for (
unsigned int i = 0; i < n_child_cell_dofs; ++i)
548 indices += n_child_cell_dofs;
555 template <
int dim,
typename Number>
556 template <
int degree>
559 const unsigned int from_level,
564 const unsigned int degree_size = (degree > -1 ? degree : fe_degree) + 1;
565 const unsigned int n_child_dofs_1d = 2 * degree_size - element_is_continuous;
566 const unsigned int n_scalar_cell_dofs =
567 Utilities::fixed_power<dim>(n_child_dofs_1d);
570 for (
unsigned int cell = 0; cell < n_owned_level_cells[from_level - 1];
573 const unsigned int n_lanes =
574 (cell + vec_size > n_owned_level_cells[from_level - 1]) ?
575 (n_owned_level_cells[from_level - 1] - cell) :
580 const unsigned int *indices =
581 &level_dof_indices[from_level][cell * n_child_cell_dofs];
582 for (
unsigned int v = 0; v < n_lanes; ++v)
584 for (
unsigned int i = 0; i < n_child_cell_dofs; ++i)
586 indices += n_child_cell_dofs;
591 degree_size * n_child_dofs_1d);
593 if (element_is_continuous)
596 degree != -1 ? 2 * degree + 1 :
599 &weights_on_refined[from_level - 1]
600 [(cell / vec_size) * three_to_dim],
603 evaluation_data.data());
604 for (
unsigned int c = 0; c < n_components; ++c)
613 prolongation_matrix_1d,
615 evaluation_data.begin() +
616 c * n_scalar_cell_dofs,
617 evaluation_data.begin() +
626 for (
unsigned int c = 0; c < n_components; ++c)
635 prolongation_matrix_1d,
637 evaluation_data.begin() +
638 c * n_scalar_cell_dofs,
639 evaluation_data.begin() +
648 for (
unsigned int v = 0; v < n_lanes; ++v)
650 const unsigned int shift =
651 internal::MGTransfer::compute_shift_within_children<dim>(
652 parent_child_connect[from_level - 1][cell + v].
second,
653 fe_degree + 1 - element_is_continuous,
656 parent_child_connect[from_level - 1][cell + v].
first *
658 n_child_cell_dofs - 1,
659 level_dof_indices[from_level - 1].size());
660 const unsigned int *indices =
661 &level_dof_indices[from_level - 1]
662 [parent_child_connect[from_level - 1][cell + v]
666 for (
unsigned int c = 0, m = 0; c < n_components; ++c)
669 for (std::vector<unsigned short>::const_iterator i =
670 dirichlet_indices[from_level - 1][cell + v].
begin();
671 i != dirichlet_indices[from_level - 1][cell + v].end();
673 evaluation_data[*i][v] = 0.;
675 for (
unsigned int k = 0; k < (dim > 2 ? degree_size : 1); ++k)
676 for (
unsigned int j = 0; j < (dim > 1 ? degree_size : 1); ++j)
677 for (
unsigned int i = 0; i < degree_size; ++i, ++m)
679 indices[c * n_scalar_cell_dofs +
680 k * n_child_dofs_1d * n_child_dofs_1d +
681 j * n_child_dofs_1d + i]) +=
682 evaluation_data[m][v];
690 template <
int dim,
typename Number>
708 template <
int dim,
typename Number>
720 template <
int dim,
typename Number>
722 const std::vector<MGConstrainedDoFs> &mg_c)
727 for (
const auto &constrained_block_dofs : mg_c)
733 template <
int dim,
typename Number>
738 Assert(this->same_for_all,
739 ExcMessage(
"This object was initialized with support for usage with "
740 "one DoFHandler for each block, but this method assumes "
741 "that the same DoFHandler is used for all the blocks!"));
744 this->matrix_free_transfer_vector[0].initialize_constraints(mg_c);
749 template <
int dim,
typename Number>
752 const std::vector<MGConstrainedDoFs> &mg_c)
754 Assert(!this->same_for_all,
755 ExcMessage(
"This object was initialized with support for using "
756 "the same DoFHandler for all the blocks, but this "
757 "method assumes that there is a separate DoFHandler "
759 AssertDimension(this->matrix_free_transfer_vector.size(), mg_c.size());
761 for (
unsigned int i = 0; i < mg_c.size(); ++i)
762 this->matrix_free_transfer_vector[i].initialize_constraints(mg_c[i]);
767 template <
int dim,
typename Number>
773 this->matrix_free_transfer_vector[0].build(dof_handler);
778 template <
int dim,
typename Number>
783 AssertDimension(this->matrix_free_transfer_vector.size(), dof_handler.size());
784 for (
unsigned int i = 0; i < dof_handler.size(); ++i)
785 this->matrix_free_transfer_vector[i].build(*dof_handler[i]);
790 template <
int dim,
typename Number>
794 std::size_t total_memory_consumption = 0;
795 for (
const auto &el : this->matrix_free_transfer_vector)
796 total_memory_consumption += el.memory_consumption();
797 return total_memory_consumption;
802 template <
int dim,
typename Number>
805 const unsigned int b)
const
807 return matrix_free_transfer_vector[
b];
812 template <
int dim,
typename Number>
816 this->matrix_free_transfer_vector.clear();
822 #include "mg_transfer_matrix_free.inst"
const Triangulation< dim, spacedim > & get_triangulation() const
bool has_level_dofs() const
void zero_out_ghost_values() const
Number local_element(const size_type local_index) const
void update_ghost_values() const
size_type locally_owned_size() const
const std::shared_ptr< const Utilities::MPI::Partitioner > & get_partitioner() const
virtual void compress(::VectorOperation::values operation) override
void reinit(const size_type size, const bool omit_zeroing_entries=false)
SmartPointer< const MGConstrainedDoFs, MGLevelGlobalTransfer< LinearAlgebra::distributed::Vector< Number > > > mg_constrained_dofs
std::vector< MGTransferMatrixFree< dim, Number > > matrix_free_transfer_vector
std::size_t memory_consumption() const
MGTransferBlockMatrixFree()=default
const MGTransferMatrixFree< dim, Number > & get_matrix_free_transfer(const unsigned int b) const override
void build(const DoFHandler< dim, dim > &dof_handler)
void initialize_constraints(const MGConstrainedDoFs &mg_constrained_dofs)
void do_restrict_add(const unsigned int from_level, LinearAlgebra::distributed::Vector< Number > &dst, const LinearAlgebra::distributed::Vector< Number > &src) const
void do_prolongate_add(const unsigned int to_level, LinearAlgebra::distributed::Vector< Number > &dst, const LinearAlgebra::distributed::Vector< Number > &src) const
virtual void prolongate(const unsigned int to_level, LinearAlgebra::distributed::Vector< Number > &dst, const LinearAlgebra::distributed::Vector< Number > &src) const override
void initialize_constraints(const MGConstrainedDoFs &mg_constrained_dofs)
void build(const DoFHandler< dim, dim > &dof_handler, const std::vector< std::shared_ptr< const Utilities::MPI::Partitioner >> &external_partitioners=std::vector< std::shared_ptr< const Utilities::MPI::Partitioner >>())
virtual void prolongate_and_add(const unsigned int to_level, LinearAlgebra::distributed::Vector< Number > &dst, const LinearAlgebra::distributed::Vector< Number > &src) const override
std::size_t memory_consumption() const
virtual void restrict_and_add(const unsigned int from_level, LinearAlgebra::distributed::Vector< Number > &dst, const LinearAlgebra::distributed::Vector< Number > &src) const override
virtual unsigned int n_global_levels() const
static constexpr std::size_t size()
#define DEAL_II_NAMESPACE_OPEN
#define DEAL_II_NAMESPACE_CLOSE
#define Assert(cond, exc)
#define AssertDimension(dim1, dim2)
#define AssertIndexRange(index, range)
static ::ExceptionBase & ExcIndexRange(std::size_t arg1, std::size_t arg2, std::size_t arg3)
static ::ExceptionBase & ExcMessage(std::string arg1)
std::enable_if_t< std::is_fundamental< T >::value, std::size_t > memory_consumption(const T &t)
SymmetricTensor< 2, dim, Number > b(const Tensor< 2, dim, Number > &F)
VectorType::value_type * begin(VectorType &V)
constexpr T pow(const T base, const int iexp)
void weight_fe_q_dofs_by_entity(const Number *weights, const unsigned int n_components, const int n_points_1d_non_template, Number *data)
std::vector< Number > prolongation_matrix_1d
unsigned int n_child_cell_dofs
bool element_is_continuous
unsigned int n_components