17 #include <deal.II/base/logstream.h> 18 #include <deal.II/base/function.h> 20 #include <deal.II/lac/la_parallel_vector.h> 21 #include <deal.II/grid/tria.h> 22 #include <deal.II/grid/tria_iterator.h> 23 #include <deal.II/dofs/dof_tools.h> 24 #include <deal.II/fe/fe.h> 25 #include <deal.II/dofs/dof_accessor.h> 26 #include <deal.II/multigrid/mg_tools.h> 27 #include <deal.II/multigrid/mg_transfer_matrix_free.h> 28 #include <deal.II/multigrid/mg_transfer_internal.h> 30 #include <deal.II/matrix_free/tensor_product_kernels.h> 34 DEAL_II_NAMESPACE_OPEN
37 template<
int dim,
typename Number>
41 element_is_continuous(false),
48 template<
int dim,
typename Number>
61 template <
int dim,
typename Number>
67 template <
int dim,
typename Number>
76 template <
int dim,
typename Number>
94 template <
int dim,
typename Number>
100 std::vector<std::vector<Number> > weights_unvectorized;
104 internal::MGTransfer::setup_transfer<dim,Number>(mg_dof,
111 weights_unvectorized,
122 for (
unsigned int i=0; i<elem_info.prolongation_matrix_1d.size(); i++)
127 const unsigned int n_levels = mg_dof.get_triangulation().n_global_levels();
129 const unsigned int n_weights_per_cell = Utilities::fixed_power<dim>(3);
131 for (
unsigned int level = 1; level<n_levels; ++level)
133 weights_on_refined[level-1].resize(((n_owned_level_cells[level-1]+vec_size-1)/vec_size)*n_weights_per_cell);
135 for (
unsigned int c=0; c<n_owned_level_cells[level-1]; ++c)
137 const unsigned int comp = c/vec_size;
138 const unsigned int v = c%vec_size;
139 for (
unsigned int i = 0; i<n_weights_per_cell; ++i)
142 weights_on_refined[level-1][comp*n_weights_per_cell+i][v] = weights_unvectorized[level-1][c*n_weights_per_cell+i];
152 template <
int dim,
typename Number>
168 this->ghosted_level_vector[to_level] = 0.;
174 do_prolongate_add<0>(to_level, this->ghosted_level_vector[to_level],
175 this->ghosted_level_vector[to_level-1]);
177 do_prolongate_add<1>(to_level, this->ghosted_level_vector[to_level],
178 this->ghosted_level_vector[to_level-1]);
180 do_prolongate_add<2>(to_level, this->ghosted_level_vector[to_level],
181 this->ghosted_level_vector[to_level-1]);
183 do_prolongate_add<3>(to_level, this->ghosted_level_vector[to_level],
184 this->ghosted_level_vector[to_level-1]);
186 do_prolongate_add<4>(to_level, this->ghosted_level_vector[to_level],
187 this->ghosted_level_vector[to_level-1]);
189 do_prolongate_add<5>(to_level, this->ghosted_level_vector[to_level],
190 this->ghosted_level_vector[to_level-1]);
192 do_prolongate_add<6>(to_level, this->ghosted_level_vector[to_level],
193 this->ghosted_level_vector[to_level-1]);
195 do_prolongate_add<7>(to_level, this->ghosted_level_vector[to_level],
196 this->ghosted_level_vector[to_level-1]);
198 do_prolongate_add<8>(to_level, this->ghosted_level_vector[to_level],
199 this->ghosted_level_vector[to_level-1]);
201 do_prolongate_add<9>(to_level, this->ghosted_level_vector[to_level],
202 this->ghosted_level_vector[to_level-1]);
204 do_prolongate_add<10>(to_level, this->ghosted_level_vector[to_level],
205 this->ghosted_level_vector[to_level-1]);
208 this->ghosted_level_vector[to_level-1]);
211 dst = this->ghosted_level_vector[to_level];
216 template <
int dim,
typename Number>
232 this->ghosted_level_vector[from_level-1] = 0.;
235 do_restrict_add<0>(from_level, this->ghosted_level_vector[from_level-1],
236 this->ghosted_level_vector[from_level]);
238 do_restrict_add<1>(from_level, this->ghosted_level_vector[from_level-1],
239 this->ghosted_level_vector[from_level]);
241 do_restrict_add<2>(from_level, this->ghosted_level_vector[from_level-1],
242 this->ghosted_level_vector[from_level]);
244 do_restrict_add<3>(from_level, this->ghosted_level_vector[from_level-1],
245 this->ghosted_level_vector[from_level]);
247 do_restrict_add<4>(from_level, this->ghosted_level_vector[from_level-1],
248 this->ghosted_level_vector[from_level]);
250 do_restrict_add<5>(from_level, this->ghosted_level_vector[from_level-1],
251 this->ghosted_level_vector[from_level]);
253 do_restrict_add<6>(from_level, this->ghosted_level_vector[from_level-1],
254 this->ghosted_level_vector[from_level]);
256 do_restrict_add<7>(from_level, this->ghosted_level_vector[from_level-1],
257 this->ghosted_level_vector[from_level]);
259 do_restrict_add<8>(from_level, this->ghosted_level_vector[from_level-1],
260 this->ghosted_level_vector[from_level]);
262 do_restrict_add<9>(from_level, this->ghosted_level_vector[from_level-1],
263 this->ghosted_level_vector[from_level]);
265 do_restrict_add<10>(from_level, this->ghosted_level_vector[from_level-1],
266 this->ghosted_level_vector[from_level]);
269 do_restrict_add<-1>(from_level, this->ghosted_level_vector[from_level-1],
270 this->ghosted_level_vector[from_level]);
273 dst += this->ghosted_level_vector[from_level-1];
280 template <
int dim,
typename Eval,
typename Number,
bool pro
longate>
282 perform_tensorized_op(
const Eval &evaluator,
298 evaluator.template values<0,prolongate,false>(t0, t2);
301 evaluator.template values<0,prolongate,false>(t0, t1);
302 evaluator.template values<1,prolongate,false>(t1, t2);
306 evaluator.template values<0,prolongate,false>(t0, t2);
307 evaluator.template values<1,prolongate,false>(t2, t1);
308 evaluator.template values<2,prolongate,false>(t1, t2);
314 t0 += Eval::dofs_per_cell;
315 t2 += Eval::n_q_points;
319 t0 += Eval::n_q_points;
320 t2 += Eval::dofs_per_cell;
325 template <
int dim,
int degree,
typename Number>
327 const unsigned int n_components,
333 const int loop_length = degree != -1 ? 2*degree+1 : 2*fe_degree+1;
334 unsigned int degree_to_3 [100];
336 for (
int i=1; i<loop_length-1; ++i)
338 degree_to_3[loop_length-1] = 2;
340 for (
int k=0; k<(dim>2 ? loop_length : 1); ++k)
341 for (
int j=0; j<(dim>1 ? loop_length : 1); ++j)
343 const unsigned int shift = 9*degree_to_3[k] + 3*degree_to_3[j];
344 data[0] *= weights[shift];
347 for (
int i=1; i<loop_length-1; ++i)
348 data[i] *= weights[shift+1];
349 data[loop_length-1] *= weights[shift+2];
357 template <
int dim,
typename Number>
358 template <
int degree>
365 const unsigned int degree_size = (degree > -1 ? degree :
fe_degree) + 1;
367 const unsigned int n_scalar_cell_dofs = Utilities::fixed_power<dim>(n_child_dofs_1d);
373 const unsigned int n_chunks = cell+vec_size > n_owned_level_cells[to_level-1] ?
374 n_owned_level_cells[to_level-1] - cell : vec_size;
377 for (
unsigned int v=0; v<n_chunks; ++v)
379 const unsigned int shift = internal::MGTransfer::compute_shift_within_children<dim>
385 for (
unsigned int k=0; k<(dim>2 ? degree_size : 1); ++k)
386 for (
unsigned int j=0; j<(dim>1 ? degree_size : 1); ++j)
387 for (
unsigned int i=0; i<degree_size; ++i, ++m)
390 k*n_child_dofs_1d*n_child_dofs_1d+
391 j*n_child_dofs_1d+i]);
400 degree_size * n_child_dofs_1d);
402 if (element_is_continuous)
410 perform_tensorized_op<dim,Evaluator,Number,true>(evaluator,
414 weight_dofs_on_child<dim,degree,Number>(&
weights_on_refined[to_level-1][(cell/vec_size)*three_to_dim],
426 perform_tensorized_op<dim,Evaluator,Number,true>(evaluator,
435 for (
unsigned int v=0; v<n_chunks; ++v)
446 template <
int dim,
typename Number>
447 template <
int degree>
454 const unsigned int degree_size = (degree > -1 ? degree :
fe_degree) + 1;
456 const unsigned int n_scalar_cell_dofs = Utilities::fixed_power<dim>(n_child_dofs_1d);
462 const unsigned int n_chunks = cell+vec_size > n_owned_level_cells[from_level-1] ?
463 n_owned_level_cells[from_level-1] - cell : vec_size;
469 for (
unsigned int v=0; v<n_chunks; ++v)
478 degree_size * n_child_dofs_1d);
480 if (element_is_continuous)
488 weight_dofs_on_child<dim,degree,Number>(&
weights_on_refined[from_level-1][(cell/vec_size)*three_to_dim],
491 perform_tensorized_op<dim,Evaluator,Number,false>(evaluator,
504 perform_tensorized_op<dim,Evaluator,Number,false>(evaluator,
511 for (
unsigned int v=0; v<n_chunks; ++v)
513 const unsigned int shift = internal::MGTransfer::compute_shift_within_children<dim>
526 for (
unsigned int k=0; k<(dim>2 ? degree_size : 1); ++k)
527 for (
unsigned int j=0; j<(dim>1 ? degree_size : 1); ++j)
528 for (
unsigned int i=0; i<degree_size; ++i, ++m)
530 k*n_child_dofs_1d*n_child_dofs_1d+
531 j*n_child_dofs_1d+i])
540 template <
int dim,
typename Number>
557 #include "mg_transfer_matrix_free.inst" 560 DEAL_II_NAMESPACE_CLOSE
static const unsigned int invalid_unsigned_int
void do_prolongate_add(const unsigned int to_level, LinearAlgebra::distributed::Vector< Number > &dst, const LinearAlgebra::distributed::Vector< Number > &src) const
#define AssertDimension(dim1, dim2)
Number local_element(const size_type local_index) const
std::vector< unsigned int > n_owned_level_cells
AlignedVector< VectorizedArray< Number > > evaluation_data
std::vector< std::vector< std::pair< unsigned int, unsigned int > > > copy_indices_global_mine
unsigned int n_components
std::vector< std::vector< std::vector< unsigned short > > > dirichlet_indices
#define AssertIndexRange(index, range)
static::ExceptionBase & ExcIndexRange(int arg1, int arg2, int arg3)
size_type local_size() const
virtual void prolongate(const unsigned int to_level, LinearAlgebra::distributed::Vector< Number > &dst, const LinearAlgebra::distributed::Vector< Number > &src) const
AlignedVector< VectorizedArray< Number > > prolongation_matrix_1d
std::vector< std::vector< std::pair< unsigned int, unsigned int > > > parent_child_connect
#define Assert(cond, exc)
std::size_t memory_consumption() const
SmartPointer< const MGConstrainedDoFs, MGLevelGlobalTransfer< LinearAlgebra::distributed::Vector< Number > > > mg_constrained_dofs
std::vector< AlignedVector< VectorizedArray< Number > > > weights_on_refined
void fill_and_communicate_copy_indices(const DoFHandler< dim, spacedim > &mg_dof)
void resize(const size_type size_in, const T &init=T())
void do_restrict_add(const unsigned int from_level, LinearAlgebra::distributed::Vector< Number > &dst, const LinearAlgebra::distributed::Vector< Number > &src) const
bool element_is_continuous
virtual ~MGTransferMatrixFree()
std::vector< std::vector< unsigned int > > level_dof_indices
std_cxx11::enable_if< std_cxx11::is_fundamental< T >::value, std::size_t >::type memory_consumption(const T &t)
void update_ghost_values() const
static::ExceptionBase & ExcNotImplemented()
unsigned int n_child_cell_dofs
virtual void restrict_and_add(const unsigned int from_level, LinearAlgebra::distributed::Vector< Number > &dst, const LinearAlgebra::distributed::Vector< Number > &src) const
unsigned int n_components(const DoFHandler< dim, spacedim > &dh)
void initialize_constraints(const MGConstrainedDoFs &mg_constrained_dofs)
MGLevelObject< LinearAlgebra::distributed::Vector< Number > > ghosted_level_vector
void build(const DoFHandler< dim, dim > &mg_dof)