1 // --------------------------------------------------------------------- 2 // 3 // Copyright (C) 2023 by the deal.II authors 4 // 5 // This file is part of the deal.II library. 6 // 7 // The deal.II library is free software; you can use it, redistribute 8 // it, and/or modify it under the terms of the GNU Lesser General 9 // Public License as published by the Free Software Foundation; either 10 // version 2.1 of the License, or (at your option) any later version. 11 // The full text of the license can be found in the file LICENSE.md at 12 // the top level directory of deal.II. 13 // 14 // Authors: Peter Munch, Martin Kronbichler 15 // 16 // --------------------------------------------------------------------- 17 18 #pragma once 19 #ifndef bps_cpu_h 20 # define bps_cpu_h 21 22 // deal.II includes 23 # include <deal.II/dofs/dof_tools.h> 24 25 # include <deal.II/fe/mapping.h> 26 27 # include <deal.II/lac/la_parallel_vector.h> 28 29 # include <deal.II/matrix_free/fe_evaluation.h> 30 # include <deal.II/matrix_free/matrix_free.h> 31 # include <deal.II/matrix_free/shape_info.h> 32 # include <deal.II/matrix_free/tools.h> 33 34 // local includes 35 # include "bps.h" 36 37 using namespace dealii; 38 39 40 41 /** 42 * Operator CPU implementation using deal.II. 43 */ 44 template <int dim, typename Number> 45 class OperatorDealii : public OperatorBase<Number, MemorySpace::Host> 46 { 47 public: 48 using VectorType = typename OperatorBase<Number, MemorySpace::Host>::VectorType; 49 50 /** 51 * Constructor. 52 */ 53 OperatorDealii(const Mapping<dim> &mapping, 54 const DoFHandler<dim> &dof_handler, 55 const AffineConstraints<Number> &constraints, 56 const Quadrature<dim> &quadrature, 57 const BPType &bp) 58 : mapping(mapping) 59 , dof_handler(dof_handler) 60 , constraints(constraints) 61 , quadrature(quadrature) 62 , bp(bp) 63 { 64 reinit(); 65 } 66 67 /** 68 * Destructor. 69 */ 70 ~OperatorDealii() = default; 71 72 /** 73 * Initialized internal data structures, particularly, MatrixFree. 74 */ 75 void 76 reinit() override 77 { 78 // configure MatrixFree 79 typename MatrixFree<dim, Number>::AdditionalData additional_data; 80 additional_data.tasks_parallel_scheme = 81 MatrixFree<dim, Number>::AdditionalData::TasksParallelScheme::none; 82 83 // create MatrixFree 84 matrix_free.reinit(mapping, dof_handler, constraints, quadrature, additional_data); 85 } 86 87 /** 88 * Matrix-vector product. 89 */ 90 void 91 vmult(VectorType &dst, const VectorType &src) const override 92 { 93 if (dof_handler.get_fe().n_components() == 1) 94 { 95 matrix_free.cell_loop(&OperatorDealii::do_cell_integral_range<1>, this, dst, src, true); 96 } 97 else 98 { 99 AssertThrow(dof_handler.get_fe().n_components() == dim, ExcInternalError()); 100 101 matrix_free.cell_loop(&OperatorDealii::do_cell_integral_range<dim>, this, dst, src, true); 102 } 103 } 104 105 /** 106 * Initialize vector. 107 */ 108 void 109 initialize_dof_vector(VectorType &vec) const override 110 { 111 matrix_free.initialize_dof_vector(vec); 112 } 113 114 /** 115 * Compute inverse of diagonal. 116 */ 117 void 118 compute_inverse_diagonal(VectorType &diagonal) const override 119 { 120 this->initialize_dof_vector(diagonal); 121 122 if (dof_handler.get_fe().n_components() == 1) 123 { 124 MatrixFreeTools::compute_diagonal(matrix_free, 125 diagonal, 126 &OperatorDealii::do_cell_integral_local<1>, 127 this); 128 } 129 else 130 { 131 AssertThrow(dof_handler.get_fe().n_components() == dim, ExcInternalError()); 132 133 MatrixFreeTools::compute_diagonal(matrix_free, 134 diagonal, 135 &OperatorDealii::do_cell_integral_local<dim>, 136 this); 137 } 138 139 for (auto &i : diagonal) 140 i = (std::abs(i) > 1.0e-10) ? (1.0 / i) : 1.0; 141 } 142 143 private: 144 /** 145 * Cell integral without vector access. 146 */ 147 template <int n_components> 148 void 149 do_cell_integral_local(FEEvaluation<dim, -1, 0, n_components, Number> &phi) const 150 { 151 if (bp <= BPType::BP2) // mass matrix 152 { 153 phi.evaluate(EvaluationFlags::values); 154 for (const auto q : phi.quadrature_point_indices()) 155 phi.submit_value(phi.get_value(q), q); 156 phi.integrate(EvaluationFlags::values); 157 } 158 else // Poisson operator 159 { 160 phi.evaluate(EvaluationFlags::gradients); 161 for (const auto q : phi.quadrature_point_indices()) 162 phi.submit_gradient(phi.get_gradient(q), q); 163 phi.integrate(EvaluationFlags::gradients); 164 } 165 } 166 167 /** 168 * Cell integral on a range of cells. 169 */ 170 template <int n_components> 171 void 172 do_cell_integral_range(const MatrixFree<dim, Number> &matrix_free, 173 VectorType &dst, 174 const VectorType &src, 175 const std::pair<unsigned int, unsigned int> &range) const 176 { 177 FEEvaluation<dim, -1, 0, n_components, Number> phi(matrix_free, range); 178 179 for (unsigned cell = range.first; cell < range.second; ++cell) 180 { 181 phi.reinit(cell); 182 phi.read_dof_values(src); // read source vector 183 do_cell_integral_local(phi); // cell integral 184 phi.distribute_local_to_global(dst); // write to destination vector 185 } 186 } 187 188 /** 189 * Mapping object passed to the constructor. 190 */ 191 const Mapping<dim> &mapping; 192 193 /** 194 * DoFHandler object passed to the constructor. 195 */ 196 const DoFHandler<dim> &dof_handler; 197 198 /** 199 * Constraints object passed to the constructor. 200 */ 201 const AffineConstraints<Number> &constraints; 202 203 /** 204 * Quadrature rule object passed to the constructor. 205 */ 206 const Quadrature<dim> &quadrature; 207 208 /** 209 * Selected BP. 210 */ 211 const BPType bp; 212 213 /** 214 * MatrixFree object. 215 */ 216 MatrixFree<dim, Number> matrix_free; 217 }; 218 219 #endif 220