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