int n_q_points_1d,
int n_components,
typename Number,
- typename VectorizedArrayType>
+ typename VectorizedArrayType,
+ typename VectorType>
void
compute_diagonal(
const MatrixFree<dim, Number, VectorizedArrayType> &matrix_free,
- LinearAlgebra::distributed::Vector<Number> & diagonal_global,
+ VectorType & diagonal_global,
const std::function<void(FEEvaluation<dim,
fe_degree,
n_q_points_1d,
int n_q_points_1d,
int n_components,
typename Number,
- typename VectorizedArrayType>
+ typename VectorizedArrayType,
+ typename VectorType>
void
compute_diagonal(
const MatrixFree<dim, Number, VectorizedArrayType> &matrix_free,
- LinearAlgebra::distributed::Vector<Number> & diagonal_global,
+ VectorType & diagonal_global,
void (CLASS::*cell_operation)(FEEvaluation<dim,
fe_degree,
n_q_points_1d,
{
this->phi.reinit(cell);
// STEP 1: get relevant information from FEEvaluation
- const unsigned int first_selected_component =
- phi.get_first_selected_component();
const auto & dof_info = phi.get_dof_info();
const unsigned int n_fe_components = dof_info.start_components.back();
const unsigned int dofs_per_component = phi.dofs_per_component;
const auto & matrix_free = phi.get_matrix_free();
+ // if we have a block vector with components with the same DoFHandler,
+ // each component is described with same set of constraints and
+ // we consider the shift in components only during access of the global
+ // vector
+ const unsigned int first_selected_component =
+ n_fe_components == 1 ? 0 : phi.get_first_selected_component();
+
const unsigned int n_lanes_filled =
matrix_free.n_active_entries_per_cell_batch(cell);
if (n_components == 1 || n_fe_components == 1)
{
- AssertDimension(n_components,
- 1); // TODO: currently no block vector supported
-
unsigned int ind_local = 0;
for (; index_indicators != next_index_indicators;
++index_indicators, ++ind_local)
}
}
}
-
// STEP 2b: sort and make unique
// sort vector
// set size locally-relevant diagonal
for (unsigned int v = 0; v < n_lanes_filled; ++v)
diagonals_local_constrained[v].assign(
- c_pools[v].row_lid_to_gid.size(), Number(0.0));
+ c_pools[v].row_lid_to_gid.size() *
+ (n_fe_components == 1 ? n_components : 1),
+ Number(0.0));
}
void
void
submit()
{
- const auto ith_column = phi.begin_dof_values();
+ // if we have a block vector with components with the same DoFHandler,
+ // we need to figure out which component and which DoF within the
+ // comonent are we currently considering
+ const unsigned int n_fe_components =
+ phi.get_dof_info().start_components.back();
+ const unsigned int comp =
+ n_fe_components == 1 ? i / phi.dofs_per_component : 0;
+ const unsigned int i_comp =
+ n_fe_components == 1 ? (i % phi.dofs_per_component) : i;
// apply local constraint matrix from left and from right:
// loop over all rows of transposed constrained matrix
const auto scale_iterator =
std::lower_bound(c_pool.col.begin() + c_pool.row[j],
c_pool.col.begin() + c_pool.row[j + 1],
- i);
+ i_comp);
// explanation: j-th row of C_e^T is empty (see above)
if (scale_iterator == c_pool.col.begin() + c_pool.row[j + 1])
continue;
// explanation: C_e^T(j,i) is zero (see above)
- if (*scale_iterator != i)
+ if (*scale_iterator != i_comp)
continue;
// apply constraint matrix from the left
Number temp = 0.0;
for (unsigned int k = c_pool.row[j]; k < c_pool.row[j + 1]; ++k)
- temp += c_pool.val[k] * ith_column[c_pool.col[k]][v];
+ temp += c_pool.val[k] *
+ phi.begin_dof_values()[comp * phi.dofs_per_component +
+ c_pool.col[k]][v];
// apply constraint matrix from the right
- diagonals_local_constrained[v][j] +=
+ diagonals_local_constrained
+ [v][j + comp * c_pools[v].row_lid_to_gid.size()] +=
temp *
c_pool.val[std::distance(c_pool.col.begin(), scale_iterator)];
}
}
}
- void
- distribute_local_to_global(
- LinearAlgebra::distributed::Vector<Number> &diagonal_global)
+ template <typename VectorType>
+ inline void
+ distribute_local_to_global(VectorType &diagonal_global)
{
// STEP 4: assembly results: add into global vector
+ const unsigned int n_fe_components =
+ phi.get_dof_info().start_components.back();
+ const unsigned int first_selected_component =
+ phi.get_first_selected_component();
+
for (unsigned int v = 0;
v < phi.get_matrix_free().n_active_entries_per_cell_batch(
phi.get_current_cell_index());
++v)
+ // if we have a block vector with components with the same
+ // DoFHandler, we need to loop over all components manully and
+ // need to apply the correct shift
for (unsigned int j = 0; j < c_pools[v].row.size() - 1; ++j)
- ::dealii::internal::vector_access_add(
- diagonal_global,
- c_pools[v].row_lid_to_gid[j],
- diagonals_local_constrained[v][j]);
+ for (unsigned int comp = 0;
+ comp < (n_fe_components == 1 ?
+ static_cast<unsigned int>(n_components) :
+ 1);
+ ++comp)
+ ::dealii::internal::vector_access_add(
+ *::dealii::internal::BlockVectorSelector<
+ VectorType,
+ IsBlockVector<VectorType>::value>::
+ get_vector_component(diagonal_global,
+ n_fe_components == 1 ?
+ comp + first_selected_component :
+ 0),
+ c_pools[v].row_lid_to_gid[j],
+ diagonals_local_constrained
+ [v][j + comp * c_pools[v].row_lid_to_gid.size()]);
}
private:
int n_q_points_1d,
int n_components,
typename Number,
- typename VectorizedArrayType>
+ typename VectorizedArrayType,
+ typename VectorType>
void
compute_diagonal(
const MatrixFree<dim, Number, VectorizedArrayType> &matrix_free,
- LinearAlgebra::distributed::Vector<Number> & diagonal_global,
+ VectorType & diagonal_global,
const std::function<void(FEEvaluation<dim,
fe_degree,
n_q_points_1d,
const unsigned int quad_no,
const unsigned int first_selected_component)
{
- using VectorType = LinearAlgebra::distributed::Vector<Number>;
-
int dummy = 0;
matrix_free.template cell_loop<VectorType, int>(
[&](const MatrixFree<dim, Number, VectorizedArrayType> &matrix_free,
- LinearAlgebra::distributed::Vector<Number> & diagonal_global,
+ VectorType & diagonal_global,
const int &,
const std::pair<unsigned int, unsigned int> &range) mutable {
FEEvaluation<dim,
int n_q_points_1d,
int n_components,
typename Number,
- typename VectorizedArrayType>
+ typename VectorizedArrayType,
+ typename VectorType>
void
compute_diagonal(
const MatrixFree<dim, Number, VectorizedArrayType> &matrix_free,
- LinearAlgebra::distributed::Vector<Number> & diagonal_global,
+ VectorType & diagonal_global,
void (CLASS::*cell_operation)(FEEvaluation<dim,
fe_degree,
n_q_points_1d,
n_q_points_1d,
n_components,
Number,
- VectorizedArrayType>(
+ VectorizedArrayType,
+ VectorType>(
matrix_free,
diagonal_global,
[&](auto &feeval) { (owning_class->*cell_operation)(feeval); },
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2019 - 2020 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE.md at
+// the top level directory of deal.II.
+//
+// ---------------------------------------------------------------------
+
+
+
+// Similar to compute_diagonal_02 but testing block vectorss.
+
+#include <deal.II/lac/la_parallel_block_vector.h>
+
+#include "compute_diagonal_util.h"
+
+using namespace dealii;
+
+template <int dim,
+ int fe_degree,
+ int n_points = fe_degree + 1,
+ int n_components = dim,
+ typename Number = double,
+ typename VectorizedArrayType = VectorizedArray<Number>>
+void
+test()
+{
+ Triangulation<dim> tria;
+ GridGenerator::hyper_ball(tria);
+ tria.refine_global(0);
+
+ const FE_Q<dim> fe_q(fe_degree);
+ const FESystem<dim> fe_system(fe_q, n_components);
+
+ DoFHandler<dim> dof_handler_q(tria);
+ dof_handler_q.distribute_dofs(fe_q);
+
+ DoFHandler<dim> dof_handler_system(tria);
+ dof_handler_system.distribute_dofs(fe_system);
+
+ AffineConstraints<Number> constraints_q;
+ DoFTools::make_hanging_node_constraints(dof_handler_q, constraints_q);
+ DoFTools::make_zero_boundary_constraints(dof_handler_q, constraints_q);
+ constraints_q.close();
+
+ constraints_q.print(std::cout);
+
+ AffineConstraints<Number> constraints_system;
+ DoFTools::make_hanging_node_constraints(dof_handler_system,
+ constraints_system);
+ DoFTools::make_zero_boundary_constraints(dof_handler_system,
+ constraints_system);
+ constraints_system.close();
+
+ constraints_system.print(std::cout);
+
+ typename MatrixFree<dim, Number, VectorizedArrayType>::AdditionalData
+ additional_data;
+ additional_data.mapping_update_flags = update_values | update_gradients;
+ additional_data.tasks_parallel_scheme =
+ MatrixFree<dim, Number, VectorizedArrayType>::AdditionalData::
+ TasksParallelScheme::none;
+
+ MappingQ<dim> mapping(1);
+ QGauss<1> quad(fe_degree + 1);
+
+ MatrixFree<dim, Number, VectorizedArrayType> matrix_free;
+ matrix_free.reinit(
+ mapping,
+ std::vector<const DoFHandler<dim> *>{&dof_handler_system, &dof_handler_q},
+ std::vector<const AffineConstraints<Number> *>{&constraints_system,
+ &constraints_q},
+ quad,
+ additional_data);
+
+ const auto kernel = [](FEEvaluation<dim,
+ fe_degree,
+ n_points,
+ n_components,
+ Number,
+ VectorizedArrayType> &phi) {
+ phi.evaluate(false, true, false);
+ for (unsigned int q = 0; q < phi.n_q_points; ++q)
+ {
+ phi.submit_symmetric_gradient(2.0 * phi.get_symmetric_gradient(q), q);
+ }
+ phi.integrate(false, true);
+ };
+
+ LinearAlgebra::distributed::Vector<Number> diagonal_global;
+ matrix_free.initialize_dof_vector(diagonal_global, 0);
+
+ MatrixFreeTools::compute_diagonal<dim,
+ fe_degree,
+ n_points,
+ n_components,
+ Number,
+ VectorizedArrayType>(matrix_free,
+ diagonal_global,
+ kernel,
+ 0);
+
+ diagonal_global.print(deallog.get_file_stream());
+
+ LinearAlgebra::distributed::BlockVector<Number> diagonal_global_block(dim);
+
+ for (unsigned int comp = 0; comp < dim; ++comp)
+ matrix_free.initialize_dof_vector(diagonal_global_block.block(comp), 1);
+
+ MatrixFreeTools::compute_diagonal<dim,
+ fe_degree,
+ n_points,
+ n_components,
+ Number,
+ VectorizedArrayType>(matrix_free,
+ diagonal_global_block,
+ kernel,
+ 1);
+
+ diagonal_global_block.print(deallog.get_file_stream());
+}
+
+int
+main(int argc, char **argv)
+{
+ Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1);
+ MPILogInitAll all;
+
+ test<2, 1>();
+}
projects / dealii.git / commitdiff