*/
namespace MatrixFreeTools
{
+ namespace internal
+ {
+ template <int dim, typename Number, bool is_face_>
+ class ComputeMatrixScratchData
+ {
+ public:
+ using FEEvalType = FEEvaluationData<dim, Number, is_face_>;
+
+ std::vector<unsigned int> dof_numbers;
+ std::vector<unsigned int> quad_numbers;
+ std::vector<unsigned int> first_selected_components;
+ std::vector<unsigned int> batch_type;
+ static const bool is_face = is_face_;
+
+ std::function<std::vector<std::unique_ptr<FEEvalType>>(
+ const std::pair<unsigned int, unsigned int> &)>
+ op_create;
+ std::function<void(std::vector<std::unique_ptr<FEEvalType>> &,
+ const unsigned int)>
+ op_reinit;
+ std::function<void(std::vector<std::unique_ptr<FEEvalType>> &)>
+ op_compute;
+ };
+ } // namespace internal
+
/**
* Modify @p additional_data so that cells are categorized
* according to their boundary IDs, making face integrals in the case of
categorize_by_boundary_ids(const Triangulation<dim> &tria,
AdditionalData &additional_data);
+
+
/**
* Compute the diagonal of a linear operator (@p diagonal_global), given
- * @p matrix_free and the local cell integral operation @p local_vmult. The
+ * @p matrix_free and the local cell integral operation @p cell_operation. The
* vector is initialized to the right size in the function.
*
* The parameters @p dof_no, @p quad_no, and @p first_selected_component are
n_q_points_1d,
n_components,
Number,
- VectorizedArrayType> &)> &local_vmult,
- const unsigned int dof_no = 0,
- const unsigned int quad_no = 0,
+ VectorizedArrayType> &)>
+ &cell_operation,
+ const unsigned int dof_no = 0,
+ const unsigned int quad_no = 0,
const unsigned int first_selected_component = 0);
+
+
/**
* Same as above but with a class and a function pointer.
*/
const unsigned int first_selected_component = 0);
+
+ /**
+ * Compute the diagonal of a linear operator (@p diagonal_global), given
+ * @p matrix_free and the local cell integral operation @p cell_operation,
+ * interior face integral operation @p face_operation, and boundary face
+ * integral operation @p boundary_operation. The
+ * vector is initialized to the right size in the function.
+ *
+ * The parameters @p dof_no, @p quad_no, and @p first_selected_component are
+ * passed to the constructor of the FEEvaluation that is internally set up.
+ */
+ template <int dim,
+ int fe_degree,
+ int n_q_points_1d,
+ int n_components,
+ typename Number,
+ typename VectorizedArrayType,
+ typename VectorType>
+ void
+ compute_diagonal(
+ const MatrixFree<dim, Number, VectorizedArrayType> &matrix_free,
+ VectorType &diagonal_global,
+ const std::function<void(FEEvaluation<dim,
+ fe_degree,
+ n_q_points_1d,
+ n_components,
+ Number,
+ VectorizedArrayType> &)>
+ &cell_operation,
+ const std::function<void(FEFaceEvaluation<dim,
+ fe_degree,
+ n_q_points_1d,
+ n_components,
+ Number,
+ VectorizedArrayType> &,
+ FEFaceEvaluation<dim,
+ fe_degree,
+ n_q_points_1d,
+ n_components,
+ Number,
+ VectorizedArrayType> &)>
+ &face_operation,
+ const std::function<void(FEFaceEvaluation<dim,
+ fe_degree,
+ n_q_points_1d,
+ n_components,
+ Number,
+ VectorizedArrayType> &)>
+ &boundary_operation,
+ const unsigned int dof_no = 0,
+ const unsigned int quad_no = 0,
+ const unsigned int first_selected_component = 0);
+
+
+
+ /**
+ * Same as above but with a class and a function pointer.
+ */
+ template <typename CLASS,
+ int dim,
+ int fe_degree,
+ int n_q_points_1d,
+ int n_components,
+ typename Number,
+ typename VectorizedArrayType,
+ typename VectorType>
+ void
+ compute_diagonal(
+ const MatrixFree<dim, Number, VectorizedArrayType> &matrix_free,
+ VectorType &diagonal_global,
+ void (CLASS::*cell_operation)(FEEvaluation<dim,
+ fe_degree,
+ n_q_points_1d,
+ n_components,
+ Number,
+ VectorizedArrayType> &) const,
+ void (CLASS::*face_operation)(FEFaceEvaluation<dim,
+ fe_degree,
+ n_q_points_1d,
+ n_components,
+ Number,
+ VectorizedArrayType> &,
+ FEFaceEvaluation<dim,
+ fe_degree,
+ n_q_points_1d,
+ n_components,
+ Number,
+ VectorizedArrayType> &)
+ const,
+ void (CLASS::*boundary_operation)(FEFaceEvaluation<dim,
+ fe_degree,
+ n_q_points_1d,
+ n_components,
+ Number,
+ VectorizedArrayType> &)
+ const,
+ const CLASS *owning_class,
+ const unsigned int dof_no = 0,
+ const unsigned int quad_no = 0,
+ const unsigned int first_selected_component = 0);
+
+
+
/**
* Compute the matrix representation of a linear operator (@p matrix), given
- * @p matrix_free and the local cell integral operation @p local_vmult.
+ * @p matrix_free and the local cell integral operation @p cell_operation.
* Constrained entries on the diagonal are set to one.
*
* The parameters @p dof_no, @p quad_no, and @p first_selected_component are
typename MatrixType>
void
compute_matrix(
- const MatrixFree<dim, Number, VectorizedArrayType> &matrix_free,
- const AffineConstraints<Number> &constraints,
- MatrixType &matrix,
+ const MatrixFree<dim, Number, VectorizedArrayType> &matrix_free,
+ const AffineConstraints<Number> &constraints,
+ MatrixType &matrix,
const std::function<void(FEEvaluation<dim,
fe_degree,
n_q_points_1d,
n_components,
Number,
- VectorizedArrayType> &)> &local_vmult,
- const unsigned int dof_no = 0,
- const unsigned int quad_no = 0,
+ VectorizedArrayType> &)>
+ &cell_operation,
+ const unsigned int dof_no = 0,
+ const unsigned int quad_no = 0,
const unsigned int first_selected_component = 0);
+
+
/**
* Same as above but with a class and a function pointer.
*/
const unsigned int first_selected_component = 0);
+ namespace internal
+ {
+ /**
+ * Compute the matrix representation of a linear operator (@p matrix), given
+ * @p matrix_free and the local cell, face and boundary integral operation.
+ */
+ template <int dim,
+ typename Number,
+ typename VectorizedArrayType,
+ typename MatrixType>
+ void
+ compute_matrix(
+ const MatrixFree<dim, Number, VectorizedArrayType> &matrix_free,
+ const AffineConstraints<Number> &constraints,
+ const internal::ComputeMatrixScratchData<dim, VectorizedArrayType, false>
+ &cell_operation,
+ const internal::ComputeMatrixScratchData<dim, VectorizedArrayType, true>
+ &face_operation,
+ const internal::ComputeMatrixScratchData<dim, VectorizedArrayType, true>
+ &boundary_operation,
+ MatrixType &matrix);
+ } // namespace internal
+
+
+
+ /**
+ * Compute the matrix representation of a linear operator (@p matrix), given
+ * @p matrix_free and the local cell integral operation @p cell_operation,
+ * interior face integral operation @p face_operation, and boundary face
+ * integal operation @p boundary_operation.
+ *
+ * The parameters @p dof_no, @p quad_no, and @p first_selected_component are
+ * passed to the constructor of the FEEvaluation that is internally set up.
+ */
+ template <int dim,
+ int fe_degree,
+ int n_q_points_1d,
+ int n_components,
+ typename Number,
+ typename VectorizedArrayType,
+ typename MatrixType>
+ void
+ compute_matrix(
+ const MatrixFree<dim, Number, VectorizedArrayType> &matrix_free,
+ const AffineConstraints<Number> &constraints,
+ MatrixType &matrix,
+ const std::function<void(FEEvaluation<dim,
+ fe_degree,
+ n_q_points_1d,
+ n_components,
+ Number,
+ VectorizedArrayType> &)>
+ &cell_operation,
+ const std::function<void(FEFaceEvaluation<dim,
+ fe_degree,
+ n_q_points_1d,
+ n_components,
+ Number,
+ VectorizedArrayType> &,
+ FEFaceEvaluation<dim,
+ fe_degree,
+ n_q_points_1d,
+ n_components,
+ Number,
+ VectorizedArrayType> &)>
+ &face_operation,
+ const std::function<void(FEFaceEvaluation<dim,
+ fe_degree,
+ n_q_points_1d,
+ n_components,
+ Number,
+ VectorizedArrayType> &)>
+ &boundary_operation,
+ const unsigned int dof_no = 0,
+ const unsigned int quad_no = 0,
+ const unsigned int first_selected_component = 0);
+
+
+
+ /**
+ * Same as above but with a class and a function pointer.
+ */
+ template <typename CLASS,
+ int dim,
+ int fe_degree,
+ int n_q_points_1d,
+ int n_components,
+ typename Number,
+ typename VectorizedArrayType,
+ typename MatrixType>
+ void
+ compute_matrix(
+ const MatrixFree<dim, Number, VectorizedArrayType> &matrix_free,
+ const AffineConstraints<Number> &constraints,
+ MatrixType &matrix,
+ void (CLASS::*cell_operation)(FEEvaluation<dim,
+ fe_degree,
+ n_q_points_1d,
+ n_components,
+ Number,
+ VectorizedArrayType> &) const,
+ void (CLASS::*face_operation)(FEFaceEvaluation<dim,
+ fe_degree,
+ n_q_points_1d,
+ n_components,
+ Number,
+ VectorizedArrayType> &,
+ FEFaceEvaluation<dim,
+ fe_degree,
+ n_q_points_1d,
+ n_components,
+ Number,
+ VectorizedArrayType> &)
+ const,
+ void (CLASS::*boundary_operation)(FEFaceEvaluation<dim,
+ fe_degree,
+ n_q_points_1d,
+ n_components,
+ Number,
+ VectorizedArrayType> &)
+ const,
+ const CLASS *owning_class,
+ const unsigned int dof_no = 0,
+ const unsigned int quad_no = 0,
+ const unsigned int first_selected_component = 0);
+
+
/**
* A wrapper around MatrixFree to help users to deal with DoFHandler
std::vector<std::pair<unsigned int, unsigned int>> inverse_lookup_origins;
};
- template <int dim,
- int fe_degree,
- int n_q_points_1d,
- int n_components,
- typename Number,
- typename VectorizedArrayType>
+ template <typename FEEvaluationType, bool is_face>
class ComputeDiagonalHelper
{
public:
- static const unsigned int n_lanes = VectorizedArrayType::size();
+ using Number = typename FEEvaluationType::number_type;
+ using VectorizedArrayType = typename FEEvaluationType::NumberType;
+
+ static const unsigned int dim = FEEvaluationType::dimension;
+ static const unsigned int n_components = FEEvaluationType::n_components;
+ static const unsigned int n_lanes = VectorizedArrayType::size();
ComputeDiagonalHelper()
: phi(nullptr)
{}
void
- initialize(FEEvaluation<dim,
- fe_degree,
- n_q_points_1d,
- n_components,
- Number,
- VectorizedArrayType> &phi)
+ initialize(FEEvaluationType &phi)
{
// if we are in hp mode and the number of unknowns changed, we must
// clear the map of entries
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);
+ this->n_lanes_filled =
+ is_face ? matrix_free.n_active_entries_per_face_batch(cell) :
+ matrix_free.n_active_entries_per_cell_batch(cell);
// STEP 2: setup CSR storage of transposed locally-relevant
// constraint matrix
std::vector<unsigned int> constraint_position;
std::vector<unsigned char> is_constrained_hn;
+ const std::array<unsigned int, n_lanes> &cells =
+ this->phi->get_cell_ids();
+
for (unsigned int v = 0; v < n_lanes_filled; ++v)
{
+ Assert(cells[v] != numbers::invalid_unsigned_int,
+ ExcInternalError());
+
const unsigned int *dof_indices;
unsigned int index_indicators, next_index_indicators;
const unsigned int start =
- (cell * n_lanes + v) * n_fe_components + first_selected_component;
+ cells[v] * n_fe_components + first_selected_component;
dof_indices =
dof_info.dof_indices.data() + dof_info.row_starts[start].first;
index_indicators = dof_info.row_starts[start].second;
[phi->get_active_fe_index()][first_selected_component])
{
const auto mask =
- dof_info.hanging_node_constraint_masks[cell * n_lanes + v];
+ dof_info.hanging_node_constraint_masks[cells[v]];
// cell has hanging nodes
if (mask != dealii::internal::MatrixFreeFunctions::
c_pools[v].row_lid_to_gid.size() *
(n_fe_components == 1 ? n_components : 1),
Number(0.0));
+
+ // check if fast path can be taken via FEEvaluation
+ bool use_fast_path = true;
+
+ for (unsigned int v = 0; v < n_lanes_filled; ++v)
+ {
+ auto &c_pool = c_pools[v];
+
+ for (unsigned int i = 0; i < c_pool.row.size() - 1; ++i)
+ {
+ if ((c_pool.row[i + 1] - c_pool.row[i]) > 1)
+ {
+ use_fast_path = false;
+ break;
+ }
+ else if (((c_pool.row[i + 1] - c_pool.row[i]) == 1) &&
+ (c_pool.val[c_pool.row[i]] != 1.0))
+ {
+ use_fast_path = false;
+ break;
+ }
+ }
+
+ if (use_fast_path == false)
+ break;
+ }
+
+ if (use_fast_path)
+ {
+ temp_values.resize(phi->dofs_per_cell);
+ return;
+ }
+ else
+ {
+ temp_values.clear();
+ }
}
void
dof_values[i] = Number(1);
}
+ void
+ zero_basis_vector()
+ {
+ VectorizedArrayType *dof_values = phi->begin_dof_values();
+ for (const unsigned int j : phi->dof_indices())
+ dof_values[j] = VectorizedArrayType();
+ }
+
void
submit()
{
+ if (!temp_values.empty())
+ {
+ temp_values[i] = phi->begin_dof_values()[i];
+ return;
+ }
+
// if we have a block vector with components with the same DoFHandler,
// we need to figure out which component and which DoF within the
// component are we currently considering
// apply local constraint matrix from left and from right:
// loop over all rows of transposed constrained matrix
- for (unsigned int v = 0;
- v < phi->get_matrix_free().n_active_entries_per_cell_batch(
- phi->get_current_cell_index());
- ++v)
+ for (unsigned int v = 0; v < n_lanes_filled; ++v)
{
const auto &c_pool = c_pools[v];
distribute_local_to_global(
std::array<VectorType *, n_components> &diagonal_global)
{
+ if (!temp_values.empty())
+ {
+ for (unsigned int j = 0; j < temp_values.size(); ++j)
+ phi->begin_dof_values()[j] = temp_values[j];
+
+ phi->distribute_local_to_global(diagonal_global);
+
+ return;
+ }
+
// STEP 4: assembly results: add into global vector
const unsigned int n_fe_components =
phi->get_dof_info().start_components.back();
- for (unsigned int v = 0;
- v < phi->get_matrix_free().n_active_entries_per_cell_batch(
- phi->get_current_cell_index());
- ++v)
+ for (unsigned int v = 0; v < n_lanes_filled; ++v)
// if we have a block vector with components with the same
// DoFHandler, we need to loop over all components manually and
// need to apply the correct shift
[v][j + comp * c_pools[v].row_lid_to_gid.size()]);
}
+ bool
+ use_fast_path() const
+ {
+ return !temp_values.empty();
+ }
+
private:
- FEEvaluation<dim,
- fe_degree,
- n_q_points_1d,
- n_components,
- Number,
- VectorizedArrayType> *phi;
+ FEEvaluationType *phi;
unsigned int dofs_per_component;
unsigned int i;
+ unsigned int n_lanes_filled;
+
std::array<internal::LocalCSR<Number>, n_lanes> c_pools;
// local storage: buffer so that we access the global vector once
locally_relevant_constraints_hn_map;
// scratch array
+ AlignedVector<VectorizedArrayType> temp_values;
AlignedVector<VectorizedArrayType> values_dofs;
};
+ template <bool is_face,
+ int dim,
+ typename Number,
+ typename VectorizedArrayType>
+ bool
+ is_fe_nothing(
+ const MatrixFree<dim, Number, VectorizedArrayType> &matrix_free,
+ const std::pair<unsigned int, unsigned int> &range,
+ const unsigned int dof_no,
+ const unsigned int quad_no,
+ const unsigned int first_selected_component,
+ const unsigned int fe_degree,
+ const unsigned int n_q_points_1d,
+ const bool is_interior_face = true)
+ {
+ const unsigned int static_n_q_points =
+ is_face ? Utilities::pow(n_q_points_1d, dim - 1) :
+ Utilities::pow(n_q_points_1d, dim);
+
+ unsigned int active_fe_index = 0;
+ if (!is_face)
+ active_fe_index = matrix_free.get_cell_active_fe_index(range);
+ else if (is_interior_face)
+ active_fe_index = matrix_free.get_face_range_category(range).first;
+ else
+ active_fe_index = matrix_free.get_face_range_category(range).second;
+
+ const auto init_data = dealii::internal::
+ extract_initialization_data<is_face, dim, Number, VectorizedArrayType>(
+ matrix_free,
+ dof_no,
+ first_selected_component,
+ quad_no,
+ fe_degree,
+ static_n_q_points,
+ active_fe_index,
+ numbers::invalid_unsigned_int /*active_quad_index*/,
+ numbers::invalid_unsigned_int /*face_type*/);
+
+ return init_data.shape_info->dofs_per_component_on_cell == 0;
+ }
+
} // namespace internal
template <int dim,
n_q_points_1d,
n_components,
Number,
- VectorizedArrayType> &)> &local_vmult,
- const unsigned int dof_no,
- const unsigned int quad_no,
+ VectorizedArrayType> &)>
+ &cell_operation,
+ const unsigned int dof_no,
+ const unsigned int quad_no,
+ const unsigned int first_selected_component)
+ {
+ compute_diagonal<dim,
+ fe_degree,
+ n_q_points_1d,
+ n_components,
+ Number,
+ VectorizedArrayType,
+ VectorType>(matrix_free,
+ diagonal_global,
+ cell_operation,
+ {},
+ {},
+ dof_no,
+ quad_no,
+ first_selected_component);
+ }
+
+ template <typename CLASS,
+ int dim,
+ int fe_degree,
+ int n_q_points_1d,
+ int n_components,
+ typename Number,
+ typename VectorizedArrayType,
+ typename VectorType>
+ void
+ compute_diagonal(
+ const MatrixFree<dim, Number, VectorizedArrayType> &matrix_free,
+ VectorType &diagonal_global,
+ void (CLASS::*cell_operation)(FEEvaluation<dim,
+ fe_degree,
+ n_q_points_1d,
+ n_components,
+ Number,
+ VectorizedArrayType> &) const,
+ const CLASS *owning_class,
+ const unsigned int dof_no,
+ const unsigned int quad_no,
+ const unsigned int first_selected_component)
+ {
+ compute_diagonal<dim,
+ fe_degree,
+ n_q_points_1d,
+ n_components,
+ Number,
+ VectorizedArrayType,
+ VectorType>(
+ matrix_free,
+ diagonal_global,
+ [&](auto &phi) { (owning_class->*cell_operation)(phi); },
+ dof_no,
+ quad_no,
+ first_selected_component);
+ }
+
+ template <int dim,
+ int fe_degree,
+ int n_q_points_1d,
+ int n_components,
+ typename Number,
+ typename VectorizedArrayType,
+ typename VectorType>
+ void
+ compute_diagonal(
+ const MatrixFree<dim, Number, VectorizedArrayType> &matrix_free,
+ VectorType &diagonal_global,
+ const std::function<void(FEEvaluation<dim,
+ fe_degree,
+ n_q_points_1d,
+ n_components,
+ Number,
+ VectorizedArrayType> &)>
+ &cell_operation,
+ const std::function<void(FEFaceEvaluation<dim,
+ fe_degree,
+ n_q_points_1d,
+ n_components,
+ Number,
+ VectorizedArrayType> &,
+ FEFaceEvaluation<dim,
+ fe_degree,
+ n_q_points_1d,
+ n_components,
+ Number,
+ VectorizedArrayType> &)>
+ &face_operation,
+ const std::function<void(FEFaceEvaluation<dim,
+ fe_degree,
+ n_q_points_1d,
+ n_components,
+ Number,
+ VectorizedArrayType> &)>
+ &boundary_operation,
+ const unsigned int dof_no,
+ const unsigned int quad_no,
const unsigned int first_selected_component)
{
int dummy = 0;
*diagonal_global_components[0], matrix_free, dof_info);
}
- using Helper = internal::ComputeDiagonalHelper<dim,
+ using Helper =
+ internal::ComputeDiagonalHelper<FEEvaluation<dim,
fe_degree,
n_q_points_1d,
n_components,
Number,
- VectorizedArrayType>;
-
- Threads::ThreadLocalStorage<Helper> scratch_data;
- matrix_free.template cell_loop<VectorType, int>(
+ VectorizedArrayType>,
+ false>;
+
+ using HelperFace =
+ internal::ComputeDiagonalHelper<FEFaceEvaluation<dim,
+ fe_degree,
+ n_q_points_1d,
+ n_components,
+ Number,
+ VectorizedArrayType>,
+ true>;
+
+ Threads::ThreadLocalStorage<Helper> scratch_data;
+ Threads::ThreadLocalStorage<HelperFace> scratch_data_m;
+ Threads::ThreadLocalStorage<HelperFace> scratch_data_p;
+
+ const auto cell_operation_wrapped =
[&](const MatrixFree<dim, Number, VectorizedArrayType> &matrix_free,
VectorType &,
const int &,
- const std::pair<unsigned int, unsigned int> &range) mutable {
+ const std::pair<unsigned int, unsigned int> &range) {
+ if (!cell_operation)
+ return;
+
+ // shortcut for FE_Nothing cells
+ if (internal::is_fe_nothing<false>(matrix_free,
+ range,
+ dof_no,
+ quad_no,
+ first_selected_component,
+ fe_degree,
+ n_q_points_1d))
+ return;
+
Helper &helper = scratch_data.get();
FEEvaluation<dim,
for (unsigned int i = 0; i < phi.dofs_per_cell; ++i)
{
helper.prepare_basis_vector(i);
- local_vmult(phi);
+ cell_operation(phi);
helper.submit();
}
helper.distribute_local_to_global(diagonal_global_components);
}
- },
- diagonal_global,
- dummy,
- false);
+ };
+
+ const auto face_operation_wrapped =
+ [&](const MatrixFree<dim, Number, VectorizedArrayType> &matrix_free,
+ VectorType &,
+ const int &,
+ const std::pair<unsigned int, unsigned int> &range) {
+ if (!face_operation)
+ return;
+
+ // shortcut for faces containing purely FE_Nothing
+ if (internal::is_fe_nothing<true>(matrix_free,
+ range,
+ dof_no,
+ quad_no,
+ first_selected_component,
+ fe_degree,
+ n_q_points_1d,
+ true) ||
+ internal::is_fe_nothing<true>(matrix_free,
+ range,
+ dof_no,
+ quad_no,
+ first_selected_component,
+ fe_degree,
+ n_q_points_1d,
+ false))
+ return;
+
+ HelperFace &helper_m = scratch_data_m.get();
+ HelperFace &helper_p = scratch_data_p.get();
+
+ FEFaceEvaluation<dim,
+ fe_degree,
+ n_q_points_1d,
+ n_components,
+ Number,
+ VectorizedArrayType>
+ phi_m(matrix_free,
+ range,
+ true,
+ dof_no,
+ quad_no,
+ first_selected_component);
+ FEFaceEvaluation<dim,
+ fe_degree,
+ n_q_points_1d,
+ n_components,
+ Number,
+ VectorizedArrayType>
+ phi_p(matrix_free,
+ range,
+ false,
+ dof_no,
+ quad_no,
+ first_selected_component);
+
+ helper_m.initialize(phi_m);
+ helper_p.initialize(phi_p);
+
+ for (unsigned int face = range.first; face < range.second; ++face)
+ {
+ helper_m.reinit(face);
+ helper_p.reinit(face);
+
+ // make check only if both adjacent cells have DoFs
+ Assert(helper_m.use_fast_path() && helper_p.use_fast_path(),
+ ExcNotImplemented());
+
+ for (unsigned int i = 0; i < phi_m.dofs_per_cell; ++i)
+ {
+ helper_m.prepare_basis_vector(i);
+ helper_p.zero_basis_vector();
+ face_operation(phi_m, phi_p);
+ helper_m.submit();
+ }
+
+ helper_m.distribute_local_to_global(diagonal_global_components);
+
+ for (unsigned int i = 0; i < phi_p.dofs_per_cell; ++i)
+ {
+ helper_m.zero_basis_vector();
+ helper_p.prepare_basis_vector(i);
+ face_operation(phi_m, phi_p);
+ helper_p.submit();
+ }
+
+ helper_p.distribute_local_to_global(diagonal_global_components);
+ }
+ };
+
+ const auto boundary_operation_wrapped =
+ [&](const MatrixFree<dim, Number, VectorizedArrayType> &matrix_free,
+ VectorType &,
+ const int &,
+ const std::pair<unsigned int, unsigned int> &range) {
+ if (!boundary_operation)
+ return;
+
+ // shortcut for faces containing purely FE_Nothing
+ if (internal::is_fe_nothing<true>(matrix_free,
+ range,
+ dof_no,
+ quad_no,
+ first_selected_component,
+ fe_degree,
+ n_q_points_1d,
+ true))
+ return;
+
+ HelperFace &helper = scratch_data_m.get();
+
+ FEFaceEvaluation<dim,
+ fe_degree,
+ n_q_points_1d,
+ n_components,
+ Number,
+ VectorizedArrayType>
+ phi(matrix_free,
+ range,
+ true,
+ dof_no,
+ quad_no,
+ first_selected_component);
+ helper.initialize(phi);
+
+ for (unsigned int face = range.first; face < range.second; ++face)
+ {
+ helper.reinit(face);
+
+ for (unsigned int i = 0; i < phi.dofs_per_cell; ++i)
+ {
+ helper.prepare_basis_vector(i);
+ boundary_operation(phi);
+ helper.submit();
+ }
+
+ helper.distribute_local_to_global(diagonal_global_components);
+ }
+ };
+
+ if (face_operation || boundary_operation)
+ matrix_free.template loop<VectorType, int>(cell_operation_wrapped,
+ face_operation_wrapped,
+ boundary_operation_wrapped,
+ diagonal_global,
+ dummy,
+ false);
+ else
+ matrix_free.template cell_loop<VectorType, int>(cell_operation_wrapped,
+ diagonal_global,
+ dummy,
+ false);
}
template <typename CLASS,
n_components,
Number,
VectorizedArrayType> &) const,
+ void (CLASS::*face_operation)(FEFaceEvaluation<dim,
+ fe_degree,
+ n_q_points_1d,
+ n_components,
+ Number,
+ VectorizedArrayType> &,
+ FEFaceEvaluation<dim,
+ fe_degree,
+ n_q_points_1d,
+ n_components,
+ Number,
+ VectorizedArrayType> &)
+ const,
+ void (CLASS::*boundary_operation)(FEFaceEvaluation<dim,
+ fe_degree,
+ n_q_points_1d,
+ n_components,
+ Number,
+ VectorizedArrayType> &)
+ const,
const CLASS *owning_class,
const unsigned int dof_no,
const unsigned int quad_no,
VectorType>(
matrix_free,
diagonal_global,
- [&](auto &feeval) { (owning_class->*cell_operation)(feeval); },
+ [&](auto &phi) { (owning_class->*cell_operation)(phi); },
+ [&](auto &phi_m, auto &phi_p) {
+ (owning_class->*face_operation)(phi_m, phi_p);
+ },
+ [&](auto &phi) { (owning_class->*boundary_operation)(phi); },
dof_no,
quad_no,
first_selected_component);
n_q_points_1d,
n_components,
Number,
- VectorizedArrayType> &)> &local_vmult,
- const unsigned int dof_no,
- const unsigned int quad_no,
+ VectorizedArrayType> &)>
+ &cell_operation,
+ const unsigned int dof_no,
+ const unsigned int quad_no,
const unsigned int first_selected_component)
{
- std::unique_ptr<AffineConstraints<typename MatrixType::value_type>>
- constraints_for_matrix;
- const AffineConstraints<typename MatrixType::value_type> &constraints =
- internal::create_new_affine_constraints_if_needed(matrix,
- constraints_in,
- constraints_for_matrix);
-
- matrix_free.template cell_loop<MatrixType, MatrixType>(
- [&](const auto &, auto &dst, const auto &, const auto range) {
- FEEvaluation<dim,
- fe_degree,
- n_q_points_1d,
- n_components,
- Number,
- VectorizedArrayType>
- integrator(
- matrix_free, range, dof_no, quad_no, first_selected_component);
-
- const unsigned int dofs_per_cell = integrator.dofs_per_cell;
+ compute_matrix<dim,
+ fe_degree,
+ n_q_points_1d,
+ n_components,
+ Number,
+ VectorizedArrayType,
+ MatrixType>(matrix_free,
+ constraints_in,
+ matrix,
+ cell_operation,
+ {},
+ {},
+ dof_no,
+ quad_no,
+ first_selected_component);
+ }
- std::vector<types::global_dof_index> dof_indices(dofs_per_cell);
- std::vector<types::global_dof_index> dof_indices_mf(dofs_per_cell);
+ template <typename CLASS,
+ int dim,
+ int fe_degree,
+ int n_q_points_1d,
+ int n_components,
+ typename Number,
+ typename VectorizedArrayType,
+ typename MatrixType>
+ void
+ compute_matrix(
+ const MatrixFree<dim, Number, VectorizedArrayType> &matrix_free,
+ const AffineConstraints<Number> &constraints,
+ MatrixType &matrix,
+ void (CLASS::*cell_operation)(FEEvaluation<dim,
+ fe_degree,
+ n_q_points_1d,
+ n_components,
+ Number,
+ VectorizedArrayType> &) const,
+ const CLASS *owning_class,
+ const unsigned int dof_no,
+ const unsigned int quad_no,
+ const unsigned int first_selected_component)
+ {
+ compute_matrix<dim,
+ fe_degree,
+ n_q_points_1d,
+ n_components,
+ Number,
+ VectorizedArrayType,
+ MatrixType>(
+ matrix_free,
+ constraints,
+ matrix,
+ [&](auto &phi) { (owning_class->*cell_operation)(phi); },
+ dof_no,
+ quad_no,
+ first_selected_component);
+ }
- std::array<FullMatrix<typename MatrixType::value_type>,
- VectorizedArrayType::size()>
- matrices;
+ template <int dim,
+ int fe_degree,
+ int n_q_points_1d,
+ int n_components,
+ typename Number,
+ typename VectorizedArrayType,
+ typename MatrixType>
+ void
+ compute_matrix(
+ const MatrixFree<dim, Number, VectorizedArrayType> &matrix_free,
+ const AffineConstraints<Number> &constraints_in,
+ MatrixType &matrix,
+ const std::function<void(FEEvaluation<dim,
+ fe_degree,
+ n_q_points_1d,
+ n_components,
+ Number,
+ VectorizedArrayType> &)>
+ &cell_operation,
+ const std::function<void(FEFaceEvaluation<dim,
+ fe_degree,
+ n_q_points_1d,
+ n_components,
+ Number,
+ VectorizedArrayType> &,
+ FEFaceEvaluation<dim,
+ fe_degree,
+ n_q_points_1d,
+ n_components,
+ Number,
+ VectorizedArrayType> &)>
+ &face_operation,
+ const std::function<void(FEFaceEvaluation<dim,
+ fe_degree,
+ n_q_points_1d,
+ n_components,
+ Number,
+ VectorizedArrayType> &)>
+ &boundary_operation,
+ const unsigned int dof_no,
+ const unsigned int quad_no,
+ const unsigned int first_selected_component)
+ {
+ using FEEvalType = FEEvaluation<dim,
+ fe_degree,
+ n_q_points_1d,
+ n_components,
+ Number,
+ VectorizedArrayType>;
+
+ using FEFaceEvalType = FEFaceEvaluation<dim,
+ fe_degree,
+ n_q_points_1d,
+ n_components,
+ Number,
+ VectorizedArrayType>;
+
+ internal::ComputeMatrixScratchData<dim, VectorizedArrayType, false>
+ data_cell;
+
+ data_cell.dof_numbers = {dof_no};
+ data_cell.quad_numbers = {dof_no};
+ data_cell.first_selected_components = {first_selected_component};
+ data_cell.batch_type = {0};
+
+ data_cell.op_create =
+ [&](const std::pair<unsigned int, unsigned int> &range) {
+ std::vector<
+ std::unique_ptr<FEEvaluationData<dim, VectorizedArrayType, false>>>
+ phi;
+
+ if (!internal::is_fe_nothing<false>(matrix_free,
+ range,
+ dof_no,
+ quad_no,
+ first_selected_component,
+ fe_degree,
+ n_q_points_1d))
+ phi.emplace_back(std::make_unique<FEEvalType>(
+ matrix_free, range, dof_no, quad_no, first_selected_component));
+
+ return phi;
+ };
- std::fill_n(matrices.begin(),
- VectorizedArrayType::size(),
- FullMatrix<typename MatrixType::value_type>(dofs_per_cell,
- dofs_per_cell));
+ data_cell.op_reinit = [](auto &phi, const unsigned batch) {
+ static_cast<FEEvalType &>(*phi[0]).reinit(batch);
+ };
- const auto lexicographic_numbering =
- matrix_free
- .get_shape_info(dof_no,
- quad_no,
- first_selected_component,
- integrator.get_active_fe_index(),
- integrator.get_active_quadrature_index())
- .lexicographic_numbering;
+ if (cell_operation)
+ data_cell.op_compute = [&](auto &phi) {
+ cell_operation(static_cast<FEEvalType &>(*phi[0]));
+ };
- for (auto cell = range.first; cell < range.second; ++cell)
+ internal::ComputeMatrixScratchData<dim, VectorizedArrayType, true>
+ data_face;
+
+ data_face.dof_numbers = {dof_no, dof_no};
+ data_face.quad_numbers = {dof_no, quad_no};
+ data_face.first_selected_components = {first_selected_component,
+ first_selected_component};
+ data_face.batch_type = {1, 2};
+
+ data_face.op_create =
+ [&](const std::pair<unsigned int, unsigned int> &range) {
+ std::vector<
+ std::unique_ptr<FEEvaluationData<dim, VectorizedArrayType, true>>>
+ phi;
+
+ if (!internal::is_fe_nothing<true>(matrix_free,
+ range,
+ dof_no,
+ quad_no,
+ first_selected_component,
+ fe_degree,
+ n_q_points_1d,
+ true) &&
+ !internal::is_fe_nothing<true>(matrix_free,
+ range,
+ dof_no,
+ quad_no,
+ first_selected_component,
+ fe_degree,
+ n_q_points_1d,
+ false))
{
- integrator.reinit(cell);
+ phi.emplace_back(
+ std::make_unique<FEFaceEvalType>(matrix_free,
+ range,
+ true,
+ dof_no,
+ quad_no,
+ first_selected_component));
+ phi.emplace_back(
+ std::make_unique<FEFaceEvalType>(matrix_free,
+ range,
+ false,
+ dof_no,
+ quad_no,
+ first_selected_component));
+ }
- const unsigned int n_filled_lanes =
- matrix_free.n_active_entries_per_cell_batch(cell);
+ return phi;
+ };
- for (unsigned int v = 0; v < n_filled_lanes; ++v)
- matrices[v] = 0.0;
+ data_face.op_reinit = [](auto &phi, const unsigned batch) {
+ static_cast<FEFaceEvalType &>(*phi[0]).reinit(batch);
+ static_cast<FEFaceEvalType &>(*phi[1]).reinit(batch);
+ };
- for (unsigned int j = 0; j < dofs_per_cell; ++j)
- {
- for (unsigned int i = 0; i < dofs_per_cell; ++i)
- integrator.begin_dof_values()[i] =
- static_cast<Number>(i == j);
+ if (face_operation)
+ data_face.op_compute = [&](auto &phi) {
+ face_operation(static_cast<FEFaceEvalType &>(*phi[0]),
+ static_cast<FEFaceEvalType &>(*phi[1]));
+ };
+
+ internal::ComputeMatrixScratchData<dim, VectorizedArrayType, true>
+ data_boundary;
+
+ data_boundary.dof_numbers = {dof_no};
+ data_boundary.quad_numbers = {dof_no};
+ data_boundary.first_selected_components = {first_selected_component};
+ data_boundary.batch_type = {1};
+
+ data_boundary
+ .op_create = [&](const std::pair<unsigned int, unsigned int> &range) {
+ std::vector<
+ std::unique_ptr<FEEvaluationData<dim, VectorizedArrayType, true>>>
+ phi;
+
+ if (!internal::is_fe_nothing<true>(matrix_free,
+ range,
+ dof_no,
+ quad_no,
+ first_selected_component,
+ fe_degree,
+ n_q_points_1d,
+ true))
+ phi.emplace_back(std::make_unique<FEFaceEvalType>(
+ matrix_free, range, true, dof_no, quad_no, first_selected_component));
+
+ return phi;
+ };
+
+ data_boundary.op_reinit = [](auto &phi, const unsigned batch) {
+ static_cast<FEFaceEvalType &>(*phi[0]).reinit(batch);
+ };
+
+ if (boundary_operation)
+ data_boundary.op_compute = [&](auto &phi) {
+ boundary_operation(static_cast<FEFaceEvalType &>(*phi[0]));
+ };
+
+ internal::compute_matrix(
+ matrix_free, constraints_in, data_cell, data_face, data_boundary, matrix);
+ }
+
+ namespace internal
+ {
+ template <int dim,
+ typename Number,
+ typename VectorizedArrayType,
+ typename MatrixType>
+ void
+ compute_matrix(
+ const MatrixFree<dim, Number, VectorizedArrayType> &matrix_free,
+ const AffineConstraints<Number> &constraints_in,
+ const internal::ComputeMatrixScratchData<dim, VectorizedArrayType, false>
+ &data_cell,
+ const internal::ComputeMatrixScratchData<dim, VectorizedArrayType, true>
+ &data_face,
+ const internal::ComputeMatrixScratchData<dim, VectorizedArrayType, true>
+ &data_boundary,
+ MatrixType &matrix)
+ {
+ std::unique_ptr<AffineConstraints<typename MatrixType::value_type>>
+ constraints_for_matrix;
+ const AffineConstraints<typename MatrixType::value_type> &constraints =
+ internal::create_new_affine_constraints_if_needed(
+ matrix, constraints_in, constraints_for_matrix);
+
+ const auto batch_operation =
+ [&matrix_free, &constraints, &matrix](
+ auto &data, const std::pair<unsigned int, unsigned int> &range) {
+ if (!data.op_compute)
+ return; // nothing to do
+
+ auto phi = data.op_create(range);
+
+ const unsigned int n_blocks = phi.size();
+
+ if (n_blocks == 0)
+ return;
+
+ Table<1, unsigned int> dofs_per_cell(n_blocks);
+
+ Table<1, std::vector<types::global_dof_index>> dof_indices(n_blocks);
+ Table<2, std::vector<types::global_dof_index>> dof_indices_mf(
+ n_blocks, VectorizedArrayType::size());
+ Table<1, std::vector<unsigned int>> lexicographic_numbering(n_blocks);
+ Table<2,
+ std::array<FullMatrix<typename MatrixType::value_type>,
+ VectorizedArrayType::size()>>
+ matrices(n_blocks, n_blocks);
+
+ for (unsigned int b = 0; b < n_blocks; ++b)
+ {
+ const auto &shape_info = matrix_free.get_shape_info(
+ data.dof_numbers[b],
+ data.quad_numbers[b],
+ data.first_selected_components[b],
+ phi[b]->get_active_fe_index(),
+ phi[b]->get_active_quadrature_index());
+
+ dofs_per_cell[b] =
+ shape_info.dofs_per_component_on_cell * shape_info.n_components;
+
+ dof_indices[b].resize(dofs_per_cell[b]);
+
+ for (unsigned int v = 0; v < VectorizedArrayType::size(); ++v)
+ dof_indices_mf[b][v].resize(dofs_per_cell[b]);
+
+ lexicographic_numbering[b] = shape_info.lexicographic_numbering;
+ }
+
+ for (unsigned int bj = 0; bj < n_blocks; ++bj)
+ for (unsigned int bi = 0; bi < n_blocks; ++bi)
+ std::fill_n(matrices[bi][bj].begin(),
+ VectorizedArrayType::size(),
+ FullMatrix<typename MatrixType::value_type>(
+ dofs_per_cell[bi], dofs_per_cell[bj]));
+
+ for (auto batch = range.first; batch < range.second; ++batch)
+ {
+ data.op_reinit(phi, batch);
+
+ const unsigned int n_filled_lanes =
+ data.is_face ?
+ matrix_free.n_active_entries_per_face_batch(batch) :
+ matrix_free.n_active_entries_per_cell_batch(batch);
+
+ for (unsigned int v = 0; v < n_filled_lanes; ++v)
+ for (unsigned int b = 0; b < n_blocks; ++b)
+ {
+ unsigned int const cell_index =
+ (data.batch_type[b] == 0) ?
+ (batch * VectorizedArrayType::size() + v) :
+ ((data.batch_type[b] == 1) ?
+ matrix_free.get_face_info(batch).cells_interior[v] :
+ matrix_free.get_face_info(batch).cells_exterior[v]);
+
+ const auto cell_iterator = matrix_free.get_cell_iterator(
+ cell_index / VectorizedArrayType::size(),
+ cell_index % VectorizedArrayType::size(),
+ data.dof_numbers[b]);
+
+ if (matrix_free.get_mg_level() !=
+ numbers::invalid_unsigned_int)
+ cell_iterator->get_mg_dof_indices(dof_indices[b]);
+ else
+ cell_iterator->get_dof_indices(dof_indices[b]);
+
+ for (unsigned int j = 0; j < dof_indices[b].size(); ++j)
+ dof_indices_mf[b][v][j] =
+ dof_indices[b][lexicographic_numbering[b][j]];
+ }
- local_vmult(integrator);
+ for (unsigned int bj = 0; bj < n_blocks; ++bj)
+ {
+ for (unsigned int j = 0; j < dofs_per_cell[bj]; ++j)
+ {
+ for (unsigned int bi = 0; bi < n_blocks; ++bi)
+ for (unsigned int i = 0; i < dofs_per_cell[bi]; ++i)
+ phi[bi]->begin_dof_values()[i] =
+ (bj == bi) ? static_cast<Number>(i == j) : 0.0;
+
+ data.op_compute(phi);
+
+ for (unsigned int bi = 0; bi < n_blocks; ++bi)
+ for (unsigned int i = 0; i < dofs_per_cell[bi]; ++i)
+ for (unsigned int v = 0; v < n_filled_lanes; ++v)
+ matrices[bi][bj][v](i, j) =
+ phi[bi]->begin_dof_values()[i][v];
+ }
- for (unsigned int i = 0; i < dofs_per_cell; ++i)
for (unsigned int v = 0; v < n_filled_lanes; ++v)
- matrices[v](i, j) = integrator.begin_dof_values()[i][v];
- }
+ for (unsigned int bi = 0; bi < n_blocks; ++bi)
+ constraints.distribute_local_to_global(
+ matrices[bi][bj][v],
+ dof_indices_mf[bi][v],
+ dof_indices_mf[bj][v],
+ matrix);
+ }
+ }
+ };
- for (unsigned int v = 0; v < n_filled_lanes; ++v)
- {
- const auto cell_v =
- matrix_free.get_cell_iterator(cell, v, dof_no);
+ const auto cell_operation_wrapped =
+ [&](const auto &, auto &, const auto &, const auto range) {
+ batch_operation(data_cell, range);
+ };
- if (matrix_free.get_mg_level() != numbers::invalid_unsigned_int)
- cell_v->get_mg_dof_indices(dof_indices);
- else
- cell_v->get_dof_indices(dof_indices);
+ const auto face_operation_wrapped =
+ [&](const auto &, auto &, const auto &, const auto range) {
+ batch_operation(data_face, range);
+ };
- for (unsigned int j = 0; j < dof_indices.size(); ++j)
- dof_indices_mf[j] = dof_indices[lexicographic_numbering[j]];
+ const auto boundary_operation_wrapped =
+ [&](const auto &, auto &, const auto &, const auto range) {
+ batch_operation(data_boundary, range);
+ };
- constraints.distribute_local_to_global(matrices[v],
- dof_indices_mf,
- dst);
- }
- }
- },
- matrix,
- matrix);
+ if (data_face.op_compute || data_boundary.op_compute)
+ {
+ matrix_free.template loop<MatrixType, MatrixType>(
+ cell_operation_wrapped,
+ face_operation_wrapped,
+ boundary_operation_wrapped,
+ matrix,
+ matrix);
+ }
+ else
+ matrix_free.template cell_loop<MatrixType, MatrixType>(
+ cell_operation_wrapped, matrix, matrix);
- matrix.compress(VectorOperation::add);
- }
+ matrix.compress(VectorOperation::add);
+ }
+ } // namespace internal
template <typename CLASS,
int dim,
n_components,
Number,
VectorizedArrayType> &) const,
+ void (CLASS::*face_operation)(FEFaceEvaluation<dim,
+ fe_degree,
+ n_q_points_1d,
+ n_components,
+ Number,
+ VectorizedArrayType> &,
+ FEFaceEvaluation<dim,
+ fe_degree,
+ n_q_points_1d,
+ n_components,
+ Number,
+ VectorizedArrayType> &)
+ const,
+ void (CLASS::*boundary_operation)(FEFaceEvaluation<dim,
+ fe_degree,
+ n_q_points_1d,
+ n_components,
+ Number,
+ VectorizedArrayType> &)
+ const,
const CLASS *owning_class,
const unsigned int dof_no,
const unsigned int quad_no,
matrix_free,
constraints,
matrix,
- [&](auto &feeval) { (owning_class->*cell_operation)(feeval); },
+ [&](auto &phi) { (owning_class->*cell_operation)(phi); },
+ [&](auto &phi_m, auto &phi_p) {
+ (owning_class->*face_operation)(phi_m, phi_p);
+ },
+ [&](auto &phi) { (owning_class->*boundary_operation)(phi); },
dof_no,
quad_no,
first_selected_component);
projects / dealii.git / commitdiff