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