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 // deal.II includes 19 #include <deal.II/dofs/dof_tools.h> 20 21 #include <deal.II/fe/mapping.h> 22 23 #include <deal.II/lac/la_parallel_vector.h> 24 25 #include <deal.II/matrix_free/fe_evaluation.h> 26 #include <deal.II/matrix_free/matrix_free.h> 27 #include <deal.II/matrix_free/tools.h> 28 29 // libCEED includes 30 #include <ceed/ceed.h> 31 32 // QFunction source 33 #include "bps-qfunctions.h" 34 35 using namespace dealii; 36 37 /** 38 * BP types. For more details, see https://ceed.exascaleproject.org/bps/. 39 */ 40 enum class BPType : unsigned int 41 { 42 BP1, 43 BP2, 44 BP3, 45 BP4, 46 BP5, 47 BP6 48 }; 49 50 51 52 /** 53 * Struct storing relevant information regarding each BP. 54 */ 55 struct BPInfo 56 { 57 BPInfo(const BPType type, const int dim, const int fe_degree) 58 : type(type) 59 , dim(dim) 60 , fe_degree(fe_degree) 61 { 62 if (type == BPType::BP1) 63 type_string = "BP1"; 64 else if (type == BPType::BP2) 65 type_string = "BP2"; 66 else if (type == BPType::BP3) 67 type_string = "BP3"; 68 else if (type == BPType::BP4) 69 type_string = "BP4"; 70 else if (type == BPType::BP5) 71 type_string = "BP5"; 72 else if (type == BPType::BP6) 73 type_string = "BP6"; 74 75 this->n_q_points_1d = (type <= BPType::BP4) ? (fe_degree + 2) : (fe_degree + 1); 76 77 this->n_components = 78 (type == BPType::BP1 || type == BPType::BP3 || type == BPType::BP5) ? 1 : dim; 79 } 80 81 82 BPType type; 83 std::string type_string; 84 unsigned int dim; 85 unsigned int fe_degree; 86 unsigned int n_q_points_1d; 87 unsigned int n_components; 88 }; 89 90 91 92 /** 93 * Base class of operators. 94 */ 95 template <typename Number> 96 class OperatorBase 97 { 98 public: 99 /** 100 * deal.II vector type 101 */ 102 using VectorType = LinearAlgebra::distributed::Vector<Number>; 103 104 /** 105 * Initialize vector. 106 */ 107 virtual void 108 reinit() = 0; 109 110 /** 111 * Perform matrix-vector product 112 */ 113 virtual void 114 vmult(VectorType &dst, const VectorType &src) const = 0; 115 116 /** 117 * Initialize vector. 118 */ 119 virtual void 120 initialize_dof_vector(VectorType &vec) const = 0; 121 122 /** 123 * Compute inverse of diagonal. 124 */ 125 virtual void 126 compute_inverse_diagonal(VectorType &diagonal) const = 0; 127 }; 128 129 130 /** 131 * Operator implementation using libCEED. 132 */ 133 template <int dim, typename Number> 134 class OperatorCeed : public OperatorBase<Number> 135 { 136 public: 137 using VectorType = typename OperatorBase<Number>::VectorType; 138 139 /** 140 * Constructor. 141 */ 142 OperatorCeed(const Mapping<dim> &mapping, 143 const DoFHandler<dim> &dof_handler, 144 const AffineConstraints<Number> &constraints, 145 const Quadrature<dim> &quadrature, 146 const BPType &bp, 147 const std::string &resource) 148 : mapping(mapping) 149 , dof_handler(dof_handler) 150 , constraints(constraints) 151 , quadrature(quadrature) 152 , bp(bp) 153 , resource(resource) 154 { 155 reinit(); 156 } 157 158 /** 159 * Destructor. 160 */ 161 ~OperatorCeed() 162 { 163 CeedOperatorDestroy(&op_apply); 164 CeedQFunctionDestroy(&qf_apply); 165 CeedQFunctionContextDestroy(&build_ctx); 166 CeedVectorDestroy(&q_data); 167 CeedElemRestrictionDestroy(&q_data_restriction); 168 CeedElemRestrictionDestroy(&sol_restriction); 169 CeedBasisDestroy(&sol_basis); 170 CeedDestroy(&ceed); 171 } 172 173 /** 174 * Initialized internal data structures, particularly, libCEED. 175 */ 176 void 177 reinit() override 178 { 179 const auto &tria = dof_handler.get_triangulation(); 180 const auto &fe = dof_handler.get_fe(); 181 182 const auto n_components = fe.n_components(); 183 184 if (bp == BPType::BP1 || bp == BPType::BP3 || bp == BPType::BP5) 185 { 186 AssertThrow(n_components == 1, ExcInternalError()); 187 } 188 else 189 { 190 AssertThrow(n_components == dim, ExcInternalError()); 191 } 192 193 // 1) create CEED instance -> "MatrixFree" 194 const char *ceed_spec = resource.c_str(); 195 CeedInit(ceed_spec, &ceed); 196 197 // 2) create shape functions -> "ShapeInfo" 198 const unsigned int fe_degree = fe.tensor_degree(); 199 const unsigned int n_q_points = quadrature.get_tensor_basis()[0].size(); 200 CeedBasisCreateTensorH1Lagrange( 201 ceed, dim, n_components, fe_degree + 1, n_q_points, CEED_GAUSS, &sol_basis); 202 203 // 3) create restriction matrix -> DoFInfo 204 unsigned int n_local_active_cells = 0; 205 206 for (const auto &cell : dof_handler.active_cell_iterators()) 207 if (cell->is_locally_owned()) 208 n_local_active_cells++; 209 210 partitioner = 211 std::make_shared<Utilities::MPI::Partitioner>(dof_handler.locally_owned_dofs(), 212 DoFTools::extract_locally_active_dofs( 213 dof_handler), 214 dof_handler.get_communicator()); 215 216 std::vector<CeedInt> indices; 217 indices.reserve(n_local_active_cells * fe.n_dofs_per_cell() / n_components); 218 219 const auto dof_mapping = FETools::lexicographic_to_hierarchic_numbering<dim>(fe_degree); 220 221 std::vector<types::global_dof_index> local_indices(fe.n_dofs_per_cell()); 222 223 for (const auto &cell : dof_handler.active_cell_iterators()) 224 if (cell->is_locally_owned()) 225 { 226 cell->get_dof_indices(local_indices); 227 228 for (const auto i : dof_mapping) 229 indices.emplace_back( 230 partitioner->global_to_local(local_indices[fe.component_to_system_index(0, i)]) / 231 n_components); 232 } 233 234 CeedElemRestrictionCreate(ceed, 235 n_local_active_cells, 236 fe.n_dofs_per_cell() / n_components, 237 n_components, 238 std::max<unsigned int>(this->extended_local_size() / n_components, 1), 239 this->extended_local_size(), 240 CEED_MEM_HOST, 241 CEED_COPY_VALUES, 242 indices.data(), 243 &sol_restriction); 244 245 // 4) create mapping -> MappingInfo 246 const unsigned int n_components_metric = (bp <= BPType::BP2) ? 1 : (dim * (dim + 1) / 2); 247 248 this->weights = compute_metric_data(ceed, mapping, tria, quadrature, bp); 249 250 strides = {{1, 251 static_cast<int>(quadrature.size()), 252 static_cast<int>(quadrature.size() * n_components_metric)}}; 253 CeedVectorCreate(ceed, weights.size(), &q_data); 254 CeedVectorSetArray(q_data, CEED_MEM_HOST, CEED_USE_POINTER, weights.data()); 255 CeedElemRestrictionCreateStrided(ceed, 256 n_local_active_cells, 257 quadrature.size(), 258 n_components_metric, 259 weights.size(), 260 strides.data(), 261 &q_data_restriction); 262 263 build_ctx_data.dim = dim; 264 build_ctx_data.space_dim = dim; 265 266 CeedQFunctionContextCreate(ceed, &build_ctx); 267 CeedQFunctionContextSetData( 268 build_ctx, CEED_MEM_HOST, CEED_USE_POINTER, sizeof(build_ctx_data), &build_ctx_data); 269 270 // 5) create q operation 271 if (bp == BPType::BP1) 272 CeedQFunctionCreateInterior(ceed, 1, f_apply_mass, f_apply_mass_loc, &qf_apply); 273 else if (bp == BPType::BP2) 274 CeedQFunctionCreateInterior(ceed, 1, f_apply_mass_vec, f_apply_mass_vec_loc, &qf_apply); 275 else if (bp == BPType::BP3 || bp == BPType::BP5) 276 CeedQFunctionCreateInterior(ceed, 1, f_apply_poisson, f_apply_poisson_loc, &qf_apply); 277 else if (bp == BPType::BP4 || bp == BPType::BP6) 278 CeedQFunctionCreateInterior(ceed, 1, f_apply_poisson_vec, f_apply_poisson_vec_loc, &qf_apply); 279 else 280 AssertThrow(false, ExcInternalError()); 281 282 if (bp <= BPType::BP2) 283 CeedQFunctionAddInput(qf_apply, "u", n_components, CEED_EVAL_INTERP); 284 else 285 CeedQFunctionAddInput(qf_apply, "u", dim * n_components, CEED_EVAL_GRAD); 286 287 CeedQFunctionAddInput(qf_apply, "qdata", n_components_metric, CEED_EVAL_NONE); 288 289 if (bp <= BPType::BP2) 290 CeedQFunctionAddOutput(qf_apply, "v", n_components, CEED_EVAL_INTERP); 291 else 292 CeedQFunctionAddOutput(qf_apply, "v", dim * n_components, CEED_EVAL_GRAD); 293 294 CeedQFunctionSetContext(qf_apply, build_ctx); 295 296 // 6) put everything together 297 CeedOperatorCreate(ceed, qf_apply, CEED_QFUNCTION_NONE, CEED_QFUNCTION_NONE, &op_apply); 298 299 CeedOperatorSetField(op_apply, "u", sol_restriction, sol_basis, CEED_VECTOR_ACTIVE); 300 CeedOperatorSetField(op_apply, "qdata", q_data_restriction, CEED_BASIS_COLLOCATED, q_data); 301 CeedOperatorSetField(op_apply, "v", sol_restriction, sol_basis, CEED_VECTOR_ACTIVE); 302 } 303 304 /** 305 * Perform matrix-vector product. 306 */ 307 void 308 vmult(VectorType &dst, const VectorType &src) const override 309 { 310 // communicate: update ghost values 311 src.update_ghost_values(); 312 313 if (dof_handler.get_fe().n_components() == 1) 314 { 315 // create libCEED view on deal.II vectors 316 VectorTypeCeed src_ceed(ceed, src); 317 VectorTypeCeed dst_ceed(ceed, dst); 318 319 // apply operator 320 CeedOperatorApply(op_apply, src_ceed(), dst_ceed(), CEED_REQUEST_IMMEDIATE); 321 } 322 else // TODO: needed for multiple components 323 { 324 // allocate space for block vectors 325 src_tmp.reinit(this->extended_local_size(), true); 326 dst_tmp.reinit(this->extended_local_size(), true); 327 328 copy_to_block_vector(src_tmp, src); // copy to block vector 329 330 // create libCEED view on deal.II vectors 331 VectorTypeCeed src_ceed(ceed, src_tmp); 332 VectorTypeCeed dst_ceed(ceed, dst_tmp); 333 334 // apply operator 335 CeedOperatorApply(op_apply, src_ceed(), dst_ceed(), CEED_REQUEST_IMMEDIATE); 336 337 dst_ceed.sync_to_host(); // pull libCEED data back to host 338 copy_from_block_vector(dst, dst_tmp); // copy from block vector 339 } 340 341 // communicate: compress 342 src.zero_out_ghost_values(); 343 dst.compress(VectorOperation::add); 344 345 // apply constraints: we assume homogeneous DBC 346 constraints.set_zero(dst); 347 } 348 349 /** 350 * Initialized vector. 351 */ 352 void 353 initialize_dof_vector(VectorType &vec) const override 354 { 355 vec.reinit(partitioner); 356 } 357 358 /** 359 * Compute inverse of diagonal. 360 */ 361 void 362 compute_inverse_diagonal(VectorType &diagonal) const override 363 { 364 this->initialize_dof_vector(diagonal); 365 366 VectorTypeCeed diagonal_ceed(ceed, diagonal); 367 368 CeedOperatorLinearAssembleDiagonal(op_apply, diagonal_ceed(), CEED_REQUEST_IMMEDIATE); 369 370 const unsigned int n_components = dof_handler.get_fe().n_components(); 371 372 if (n_components > 1) // TODO: needed for multiple components 373 { 374 VectorType tmp(diagonal); 375 376 copy_from_block_vector(tmp, diagonal); 377 378 std::swap(tmp, diagonal); 379 } 380 381 diagonal.compress(VectorOperation::add); 382 383 for (auto &i : diagonal) 384 i = (std::abs(i) > 1.0e-10) ? (1.0 / i) : 1.0; 385 } 386 387 private: 388 /** 389 * Wrapper around a deal.II vector to create a libCEED vector view. 390 */ 391 class VectorTypeCeed 392 { 393 public: 394 /** 395 * Constructor. 396 */ 397 VectorTypeCeed(const Ceed &ceed, const VectorType &vec) 398 { 399 const unsigned int n_dofs = 400 vec.get_partitioner()->locally_owned_size() + vec.get_partitioner()->n_ghost_indices(); 401 402 CeedVectorCreate(ceed, n_dofs, &vec_ceed); 403 CeedVectorSetArray(vec_ceed, CEED_MEM_HOST, CEED_USE_POINTER, vec.get_values()); 404 } 405 406 /** 407 * Return libCEED vector view. 408 */ 409 CeedVector & 410 operator()() 411 { 412 return vec_ceed; 413 } 414 415 /** 416 * Sync memory from device to host. 417 */ 418 void 419 sync_to_host() 420 { 421 CeedVectorSyncArray(vec_ceed, CEED_MEM_HOST); 422 } 423 424 /** 425 * Destructor: destroy vector view. 426 */ 427 ~VectorTypeCeed() 428 { 429 CeedScalar *ptr; 430 CeedVectorTakeArray(vec_ceed, CEED_MEM_HOST, &ptr); 431 CeedVectorDestroy(&vec_ceed); 432 } 433 434 private: 435 /** 436 * libCEED vector view. 437 */ 438 CeedVector vec_ceed; 439 }; 440 441 /** 442 * Copy from block vector. 443 * 444 * @note Only needed for multiple components. 445 */ 446 void 447 copy_from_block_vector(VectorType &dst, const VectorType &src) const 448 { 449 const unsigned int scalar_size = this->extended_local_size() / dim; 450 451 for (unsigned int i = 0; i < scalar_size; ++i) 452 for (unsigned int j = 0; j < dim; ++j) 453 dst.get_values()[j + i * dim] = src.get_values()[j * scalar_size + i]; 454 } 455 456 /** 457 * Copy to block vector. 458 * 459 * @note Only needed for multiple components. 460 */ 461 void 462 copy_to_block_vector(VectorType &dst, const VectorType &src) const 463 { 464 const unsigned int scalar_size = this->extended_local_size() / dim; 465 466 for (unsigned int i = 0; i < scalar_size; ++i) 467 for (unsigned int j = 0; j < dim; ++j) 468 dst.get_values()[j * scalar_size + i] = src.get_values()[j + i * dim]; 469 } 470 471 /** 472 * Number of locally active DoFs. 473 */ 474 unsigned int 475 extended_local_size() const 476 { 477 return partitioner->locally_owned_size() + partitioner->n_ghost_indices(); 478 } 479 480 /** 481 * Compute metric data: Jacobian, ... 482 */ 483 static std::vector<double> 484 compute_metric_data(const Ceed &ceed, 485 const Mapping<dim> &mapping, 486 const Triangulation<dim> &tria, 487 const Quadrature<dim> &quadrature, 488 const BPType bp) 489 { 490 std::vector<double> weights; 491 492 if (false) 493 { 494 FE_Nothing<dim> dummy_fe; 495 FEValues<dim> fe_values(mapping, dummy_fe, quadrature, update_JxW_values); 496 497 for (const auto &cell : tria.active_cell_iterators()) 498 if (cell->is_locally_owned()) 499 { 500 fe_values.reinit(cell); 501 502 for (const auto q : fe_values.quadrature_point_indices()) 503 weights.emplace_back(fe_values.JxW(q)); 504 } 505 506 return weights; 507 } 508 509 CeedBasis geo_basis; 510 CeedVector q_data; 511 CeedElemRestriction q_data_restriction; 512 CeedVector node_coords; 513 CeedElemRestriction geo_restriction; 514 CeedQFunctionContext build_ctx; 515 CeedQFunction qf_build; 516 CeedOperator op_build; 517 518 const unsigned int n_q_points = quadrature.get_tensor_basis()[0].size(); 519 520 const unsigned int n_components = (bp <= BPType::BP2) ? 1 : (dim * (dim + 1) / 2); 521 522 const auto mapping_q = dynamic_cast<const MappingQ<dim> *>(&mapping); 523 524 AssertThrow(mapping_q, ExcMessage("Wrong mapping!")); 525 526 const unsigned int fe_degree = mapping_q->get_degree(); 527 528 CeedBasisCreateTensorH1Lagrange( 529 ceed, dim, dim, fe_degree + 1, n_q_points, CEED_GAUSS, &geo_basis); 530 531 unsigned int n_local_active_cells = 0; 532 533 for (const auto &cell : tria.active_cell_iterators()) 534 if (cell->is_locally_owned()) 535 n_local_active_cells++; 536 537 std::vector<double> geo_support_points; 538 std::vector<CeedInt> geo_indices; 539 540 FE_Q<dim> geo_fe(fe_degree); 541 542 DoFHandler<dim> geo_dof_handler(tria); 543 geo_dof_handler.distribute_dofs(geo_fe); 544 545 const auto geo_partitioner = 546 std::make_shared<Utilities::MPI::Partitioner>(geo_dof_handler.locally_owned_dofs(), 547 DoFTools::extract_locally_active_dofs( 548 geo_dof_handler), 549 geo_dof_handler.get_communicator()); 550 551 geo_indices.reserve(n_local_active_cells * geo_fe.n_dofs_per_cell()); 552 553 const auto dof_mapping = FETools::lexicographic_to_hierarchic_numbering<dim>(fe_degree); 554 555 FEValues<dim> fe_values(mapping, 556 geo_fe, 557 geo_fe.get_unit_support_points(), 558 update_quadrature_points); 559 560 std::vector<types::global_dof_index> local_indices(geo_fe.n_dofs_per_cell()); 561 562 const unsigned int n_points = 563 geo_partitioner->locally_owned_size() + geo_partitioner->n_ghost_indices(); 564 565 geo_support_points.resize(dim * n_points); 566 567 for (const auto &cell : geo_dof_handler.active_cell_iterators()) 568 if (cell->is_locally_owned()) 569 { 570 fe_values.reinit(cell); 571 cell->get_dof_indices(local_indices); 572 573 for (const auto i : dof_mapping) 574 { 575 const auto index = geo_partitioner->global_to_local(local_indices[i]); 576 geo_indices.emplace_back(index); 577 578 const auto point = fe_values.quadrature_point(i); 579 580 for (unsigned int d = 0; d < dim; ++d) 581 geo_support_points[index + d * n_points] = point[d]; 582 } 583 } 584 585 weights.resize(n_local_active_cells * quadrature.size() * n_components); 586 587 CeedInt strides[3] = {1, 588 static_cast<int>(quadrature.size()), 589 static_cast<int>(quadrature.size() * n_components)}; 590 591 CeedVectorCreate(ceed, weights.size(), &q_data); 592 CeedVectorSetArray(q_data, CEED_MEM_HOST, CEED_USE_POINTER, weights.data()); 593 CeedElemRestrictionCreateStrided(ceed, 594 n_local_active_cells, 595 quadrature.size(), 596 n_components, 597 weights.size(), 598 strides, 599 &q_data_restriction); 600 601 CeedVectorCreate(ceed, geo_support_points.size(), &node_coords); 602 CeedVectorSetArray(node_coords, CEED_MEM_HOST, CEED_USE_POINTER, geo_support_points.data()); 603 604 CeedElemRestrictionCreate(ceed, 605 n_local_active_cells, 606 geo_fe.n_dofs_per_cell(), 607 dim, 608 std::max<unsigned int>(geo_support_points.size() / dim, 1), 609 geo_support_points.size(), 610 CEED_MEM_HOST, 611 CEED_COPY_VALUES, 612 geo_indices.data(), 613 &geo_restriction); 614 615 BuildContext build_ctx_data; 616 build_ctx_data.dim = dim; 617 build_ctx_data.space_dim = dim; 618 619 CeedQFunctionContextCreate(ceed, &build_ctx); 620 CeedQFunctionContextSetData( 621 build_ctx, CEED_MEM_HOST, CEED_USE_POINTER, sizeof(build_ctx_data), &build_ctx_data); 622 623 // 5) create q operation 624 if (bp <= BPType::BP2) 625 CeedQFunctionCreateInterior(ceed, 1, f_build_mass, f_build_mass_loc, &qf_build); 626 else 627 CeedQFunctionCreateInterior(ceed, 1, f_build_poisson, f_build_poisson_loc, &qf_build); 628 629 CeedQFunctionAddInput(qf_build, "geo", dim * dim, CEED_EVAL_GRAD); 630 CeedQFunctionAddInput(qf_build, "weights", 1, CEED_EVAL_WEIGHT); 631 CeedQFunctionAddOutput(qf_build, "qdata", n_components, CEED_EVAL_NONE); 632 CeedQFunctionSetContext(qf_build, build_ctx); 633 634 // 6) put everything together 635 CeedOperatorCreate(ceed, qf_build, CEED_QFUNCTION_NONE, CEED_QFUNCTION_NONE, &op_build); 636 CeedOperatorSetField(op_build, "geo", geo_restriction, geo_basis, CEED_VECTOR_ACTIVE); 637 CeedOperatorSetField( 638 op_build, "weights", CEED_ELEMRESTRICTION_NONE, geo_basis, CEED_VECTOR_NONE); 639 CeedOperatorSetField( 640 op_build, "qdata", q_data_restriction, CEED_BASIS_COLLOCATED, CEED_VECTOR_ACTIVE); 641 642 CeedOperatorApply(op_build, node_coords, q_data, CEED_REQUEST_IMMEDIATE); 643 644 CeedOperatorDestroy(&op_build); 645 CeedQFunctionDestroy(&qf_build); 646 CeedQFunctionContextDestroy(&build_ctx); 647 CeedElemRestrictionDestroy(&geo_restriction); 648 CeedVectorDestroy(&node_coords); 649 CeedElemRestrictionDestroy(&q_data_restriction); 650 CeedVectorSyncArray(q_data, CEED_MEM_HOST); 651 CeedVectorDestroy(&q_data); 652 CeedBasisDestroy(&geo_basis); 653 654 return weights; 655 } 656 657 /** 658 * Mapping object passed to the constructor. 659 */ 660 const Mapping<dim> &mapping; 661 662 /** 663 * DoFHandler object passed to the constructor. 664 */ 665 const DoFHandler<dim> &dof_handler; 666 667 /** 668 * Constraints object passed to the constructor. 669 */ 670 const AffineConstraints<Number> &constraints; 671 672 /** 673 * Quadrature rule object passed to the constructor. 674 */ 675 const Quadrature<dim> &quadrature; 676 677 /** 678 * Selected BP. 679 */ 680 const BPType bp; 681 682 /** 683 * Resource name. 684 */ 685 const std::string resource; 686 687 /** 688 * Partitioner for distributed vectors. 689 */ 690 std::shared_ptr<Utilities::MPI::Partitioner> partitioner; 691 692 /** 693 * libCEED data structures. 694 */ 695 Ceed ceed; 696 CeedBasis sol_basis; 697 CeedElemRestriction sol_restriction; 698 CeedElemRestriction q_data_restriction; 699 std::vector<double> weights; 700 CeedVector q_data; 701 std::array<CeedInt, 3> strides; 702 BuildContext build_ctx_data; 703 CeedQFunctionContext build_ctx; 704 CeedQFunction qf_apply; 705 CeedOperator op_apply; 706 707 /** 708 * Temporal (tempral) vectors. 709 * 710 * @note Only needed for multiple components. 711 */ 712 mutable VectorType src_tmp; 713 mutable VectorType dst_tmp; 714 }; 715 716 717 718 template <int dim, typename Number> 719 class OperatorDealii : public OperatorBase<Number> 720 { 721 public: 722 using VectorType = typename OperatorBase<Number>::VectorType; 723 724 /** 725 * Constructor. 726 */ 727 OperatorDealii(const Mapping<dim> &mapping, 728 const DoFHandler<dim> &dof_handler, 729 const AffineConstraints<Number> &constraints, 730 const Quadrature<dim> &quadrature, 731 const BPType &bp) 732 : mapping(mapping) 733 , dof_handler(dof_handler) 734 , constraints(constraints) 735 , quadrature(quadrature) 736 , bp(bp) 737 { 738 reinit(); 739 } 740 741 /** 742 * Destructor. 743 */ 744 ~OperatorDealii() = default; 745 746 /** 747 * Initialized internal data structures, particularly, MatrixFree. 748 */ 749 void 750 reinit() override 751 { 752 // configure MatrixFree 753 typename MatrixFree<dim, Number>::AdditionalData additional_data; 754 additional_data.tasks_parallel_scheme = 755 MatrixFree<dim, Number>::AdditionalData::TasksParallelScheme::none; 756 757 // create MatrixFree 758 matrix_free.reinit(mapping, dof_handler, constraints, quadrature, additional_data); 759 } 760 761 /** 762 * Matrix-vector product. 763 */ 764 void 765 vmult(VectorType &dst, const VectorType &src) const override 766 { 767 if (dof_handler.get_fe().n_components() == 1) 768 { 769 matrix_free.cell_loop(&OperatorDealii::do_cell_integral_range<1>, this, dst, src, true); 770 } 771 else 772 { 773 AssertThrow(dof_handler.get_fe().n_components() == dim, ExcInternalError()); 774 775 matrix_free.cell_loop(&OperatorDealii::do_cell_integral_range<dim>, this, dst, src, true); 776 } 777 } 778 779 /** 780 * Initialize vector. 781 */ 782 void 783 initialize_dof_vector(VectorType &vec) const override 784 { 785 matrix_free.initialize_dof_vector(vec); 786 } 787 788 /** 789 * Compute inverse of diagonal. 790 */ 791 void 792 compute_inverse_diagonal(VectorType &diagonal) const override 793 { 794 this->initialize_dof_vector(diagonal); 795 796 if (dof_handler.get_fe().n_components() == 1) 797 { 798 MatrixFreeTools::compute_diagonal(matrix_free, 799 diagonal, 800 &OperatorDealii::do_cell_integral_local<1>, 801 this); 802 } 803 else 804 { 805 AssertThrow(dof_handler.get_fe().n_components() == dim, ExcInternalError()); 806 807 MatrixFreeTools::compute_diagonal(matrix_free, 808 diagonal, 809 &OperatorDealii::do_cell_integral_local<dim>, 810 this); 811 } 812 813 for (auto &i : diagonal) 814 i = (std::abs(i) > 1.0e-10) ? (1.0 / i) : 1.0; 815 } 816 817 private: 818 /** 819 * Cell integral without vector access. 820 */ 821 template <int n_components> 822 void 823 do_cell_integral_local(FEEvaluation<dim, -1, 0, n_components, Number> &phi) const 824 { 825 if (bp <= BPType::BP2) // mass matrix 826 { 827 phi.evaluate(EvaluationFlags::values); 828 for (const auto q : phi.quadrature_point_indices()) 829 phi.submit_value(phi.get_value(q), q); 830 phi.integrate(EvaluationFlags::values); 831 } 832 else // Poisson operator 833 { 834 phi.evaluate(EvaluationFlags::gradients); 835 for (const auto q : phi.quadrature_point_indices()) 836 phi.submit_gradient(phi.get_gradient(q), q); 837 phi.integrate(EvaluationFlags::gradients); 838 } 839 } 840 841 /** 842 * Cell integral on a range of cells. 843 */ 844 template <int n_components> 845 void 846 do_cell_integral_range(const MatrixFree<dim, Number> &matrix_free, 847 VectorType &dst, 848 const VectorType &src, 849 const std::pair<unsigned int, unsigned int> &range) const 850 { 851 FEEvaluation<dim, -1, 0, n_components, Number> phi(matrix_free, range); 852 853 for (unsigned cell = range.first; cell < range.second; ++cell) 854 { 855 phi.reinit(cell); 856 phi.read_dof_values(src); // read source vector 857 do_cell_integral_local(phi); // cell integral 858 phi.distribute_local_to_global(dst); // write to destination vector 859 } 860 } 861 862 /** 863 * Mapping object passed to the constructor. 864 */ 865 const Mapping<dim> &mapping; 866 867 /** 868 * DoFHandler object passed to the constructor. 869 */ 870 const DoFHandler<dim> &dof_handler; 871 872 /** 873 * Constraints object passed to the constructor. 874 */ 875 const AffineConstraints<Number> &constraints; 876 877 /** 878 * Quadrature rule object passed to the constructor. 879 */ 880 const Quadrature<dim> &quadrature; 881 882 /** 883 * Selected BP. 884 */ 885 const BPType bp; 886 887 /** 888 * MatrixFree object. 889 */ 890 MatrixFree<dim, Number> matrix_free; 891 }; 892