1 // --------------------------------------------------------------------- 2 // 3 // Copyright (C) 2016 - 2020 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 // --------------------------------------------------------------------- 15 16 #ifndef dealii_matrix_creator_templates_h 17 #define dealii_matrix_creator_templates_h 18 19 #include <deal.II/base/config.h> 20 21 #include <deal.II/base/function.h> 22 #include <deal.II/base/geometry_info.h> 23 #include <deal.II/base/quadrature.h> 24 #include <deal.II/base/work_stream.h> 25 26 #include <deal.II/dofs/dof_accessor.h> 27 #include <deal.II/dofs/dof_handler.h> 28 #include <deal.II/dofs/dof_tools.h> 29 30 #include <deal.II/fe/fe.h> 31 #include <deal.II/fe/fe_values.h> 32 #include <deal.II/fe/mapping_q1.h> 33 34 #include <deal.II/grid/tria_iterator.h> 35 36 #include <deal.II/hp/fe_values.h> 37 #include <deal.II/hp/mapping_collection.h> 38 39 #include <deal.II/lac/block_sparse_matrix.h> 40 #include <deal.II/lac/block_vector.h> 41 #include <deal.II/lac/full_matrix.h> 42 #include <deal.II/lac/sparse_matrix.h> 43 #include <deal.II/lac/vector.h> 44 45 #include <deal.II/numerics/matrix_tools.h> 46 47 #ifdef DEAL_II_WITH_PETSC 48 # include <deal.II/lac/petsc_block_sparse_matrix.h> 49 # include <deal.II/lac/petsc_sparse_matrix.h> 50 # include <deal.II/lac/petsc_vector.h> 51 #endif 52 53 #ifdef DEAL_II_WITH_TRILINOS 54 # include <deal.II/lac/trilinos_block_sparse_matrix.h> 55 # include <deal.II/lac/trilinos_parallel_block_vector.h> 56 # include <deal.II/lac/trilinos_sparse_matrix.h> 57 # include <deal.II/lac/trilinos_vector.h> 58 #endif 59 60 61 #include <algorithm> 62 #include <cmath> 63 #include <set> 64 65 66 DEAL_II_NAMESPACE_OPEN 67 namespace MatrixCreator 68 { 69 namespace internal 70 { 71 namespace AssemblerData 72 { 73 template <int dim, int spacedim, typename number> 74 struct Scratch 75 { ScratchScratch76 Scratch(const ::dealii::hp::FECollection<dim, spacedim> &fe, 77 const UpdateFlags update_flags, 78 const Function<spacedim, number> * coefficient, 79 const Function<spacedim, number> * rhs_function, 80 const ::dealii::hp::QCollection<dim> & quadrature, 81 const ::dealii::hp::MappingCollection<dim, spacedim> &mapping) 82 : fe_collection(fe) 83 , quadrature_collection(quadrature) 84 , mapping_collection(mapping) 85 , x_fe_values(mapping_collection, 86 fe_collection, 87 quadrature_collection, 88 update_flags) 89 , coefficient_values(quadrature_collection.max_n_quadrature_points()) 90 , coefficient_vector_values( 91 quadrature_collection.max_n_quadrature_points(), 92 dealii::Vector<number>(fe_collection.n_components())) 93 , rhs_values(quadrature_collection.max_n_quadrature_points()) 94 , rhs_vector_values(quadrature_collection.max_n_quadrature_points(), 95 dealii::Vector<number>( 96 fe_collection.n_components())) 97 , coefficient(coefficient) 98 , rhs_function(rhs_function) 99 , update_flags(update_flags) 100 {} 101 ScratchScratch102 Scratch(const Scratch &data) 103 : fe_collection(data.fe_collection) 104 , quadrature_collection(data.quadrature_collection) 105 , mapping_collection(data.mapping_collection) 106 , x_fe_values(mapping_collection, 107 fe_collection, 108 quadrature_collection, 109 data.update_flags) 110 , coefficient_values(data.coefficient_values) 111 , coefficient_vector_values(data.coefficient_vector_values) 112 , rhs_values(data.rhs_values) 113 , rhs_vector_values(data.rhs_vector_values) 114 , coefficient(data.coefficient) 115 , rhs_function(data.rhs_function) 116 , update_flags(data.update_flags) 117 {} 118 119 Scratch & 120 operator=(const Scratch &) 121 { 122 Assert(false, ExcNotImplemented()); 123 return *this; 124 } 125 126 127 const ::dealii::hp::FECollection<dim, spacedim> &fe_collection; 128 const ::dealii::hp::QCollection<dim> & quadrature_collection; 129 const ::dealii::hp::MappingCollection<dim, spacedim> 130 &mapping_collection; 131 132 ::dealii::hp::FEValues<dim, spacedim> x_fe_values; 133 134 std::vector<number> coefficient_values; 135 std::vector<dealii::Vector<number>> coefficient_vector_values; 136 std::vector<number> rhs_values; 137 std::vector<dealii::Vector<number>> rhs_vector_values; 138 139 const Function<spacedim, number> *coefficient; 140 const Function<spacedim, number> *rhs_function; 141 142 const UpdateFlags update_flags; 143 }; 144 145 146 template <typename number> 147 struct CopyData 148 { 149 std::vector<types::global_dof_index> dof_indices; 150 FullMatrix<number> cell_matrix; 151 dealii::Vector<number> cell_rhs; 152 const AffineConstraints<number> * constraints; 153 }; 154 } // namespace AssemblerData 155 156 157 template <int dim, int spacedim, typename CellIterator, typename number> 158 void mass_assembler(const CellIterator & cell,MatrixCreator::internal::AssemblerData::Scratch<dim,spacedim,number> & data,MatrixCreator::internal::AssemblerData::CopyData<number> & copy_data)159 mass_assembler( 160 const CellIterator &cell, 161 MatrixCreator::internal::AssemblerData::Scratch<dim, spacedim, number> 162 & data, 163 MatrixCreator::internal::AssemblerData::CopyData<number> ©_data) 164 { 165 data.x_fe_values.reinit(cell); 166 const FEValues<dim, spacedim> &fe_values = 167 data.x_fe_values.get_present_fe_values(); 168 169 const unsigned int dofs_per_cell = fe_values.dofs_per_cell, 170 n_q_points = fe_values.n_quadrature_points; 171 const FiniteElement<dim, spacedim> &fe = fe_values.get_fe(); 172 const unsigned int n_components = fe.n_components(); 173 174 Assert(data.rhs_function == nullptr || 175 data.rhs_function->n_components == 1 || 176 data.rhs_function->n_components == n_components, 177 ::dealii::MatrixCreator::ExcComponentMismatch()); 178 Assert(data.coefficient == nullptr || 179 data.coefficient->n_components == 1 || 180 data.coefficient->n_components == n_components, 181 ::dealii::MatrixCreator::ExcComponentMismatch()); 182 183 copy_data.cell_matrix.reinit(dofs_per_cell, dofs_per_cell); 184 copy_data.cell_rhs.reinit(dofs_per_cell); 185 186 copy_data.dof_indices.resize(dofs_per_cell); 187 cell->get_dof_indices(copy_data.dof_indices); 188 189 const bool use_rhs_function = data.rhs_function != nullptr; 190 if (use_rhs_function) 191 { 192 if (data.rhs_function->n_components == 1) 193 { 194 data.rhs_values.resize(n_q_points); 195 data.rhs_function->value_list(fe_values.get_quadrature_points(), 196 data.rhs_values); 197 } 198 else 199 { 200 data.rhs_vector_values.resize( 201 n_q_points, dealii::Vector<number>(n_components)); 202 data.rhs_function->vector_value_list( 203 fe_values.get_quadrature_points(), data.rhs_vector_values); 204 } 205 } 206 207 const bool use_coefficient = data.coefficient != nullptr; 208 if (use_coefficient) 209 { 210 if (data.coefficient->n_components == 1) 211 { 212 data.coefficient_values.resize(n_q_points); 213 data.coefficient->value_list(fe_values.get_quadrature_points(), 214 data.coefficient_values); 215 } 216 else 217 { 218 data.coefficient_vector_values.resize( 219 n_q_points, dealii::Vector<number>(n_components)); 220 data.coefficient->vector_value_list( 221 fe_values.get_quadrature_points(), 222 data.coefficient_vector_values); 223 } 224 } 225 226 227 const std::vector<double> &JxW = fe_values.get_JxW_values(); 228 for (unsigned int i = 0; i < dofs_per_cell; ++i) 229 if (fe.is_primitive()) 230 { 231 const unsigned int component_i = 232 fe.system_to_component_index(i).first; 233 const double *phi_i = &fe_values.shape_value(i, 0); 234 235 // use symmetry in the mass matrix here: 236 // just need to calculate the diagonal 237 // and half of the elements above the 238 // diagonal 239 for (unsigned int j = i; j < dofs_per_cell; ++j) 240 if ((n_components == 1) || 241 (fe.system_to_component_index(j).first == component_i)) 242 { 243 const double *phi_j = &fe_values.shape_value(j, 0); 244 number add_data = 0; 245 if (use_coefficient) 246 { 247 if (data.coefficient->n_components == 1) 248 for (unsigned int point = 0; point < n_q_points; 249 ++point) 250 add_data += 251 (data.coefficient_values[point] * phi_i[point] * 252 phi_j[point] * JxW[point]); 253 else 254 for (unsigned int point = 0; point < n_q_points; 255 ++point) 256 add_data += 257 (data.coefficient_vector_values[point]( 258 component_i) * 259 phi_i[point] * phi_j[point] * JxW[point]); 260 } 261 else 262 for (unsigned int point = 0; point < n_q_points; ++point) 263 add_data += phi_i[point] * phi_j[point] * JxW[point]; 264 265 // this is even ok for i==j, since then 266 // we just write the same value twice. 267 copy_data.cell_matrix(i, j) = add_data; 268 copy_data.cell_matrix(j, i) = add_data; 269 } 270 271 if (use_rhs_function) 272 { 273 number add_data = 0; 274 if (data.rhs_function->n_components == 1) 275 for (unsigned int point = 0; point < n_q_points; ++point) 276 add_data += 277 data.rhs_values[point] * phi_i[point] * JxW[point]; 278 else 279 for (unsigned int point = 0; point < n_q_points; ++point) 280 add_data += data.rhs_vector_values[point](component_i) * 281 phi_i[point] * JxW[point]; 282 copy_data.cell_rhs(i) = add_data; 283 } 284 } 285 else 286 { 287 // non-primitive vector-valued FE, using 288 // symmetry again 289 for (unsigned int j = i; j < dofs_per_cell; ++j) 290 { 291 number add_data = 0; 292 for (unsigned int comp_i = 0; comp_i < n_components; ++comp_i) 293 if (fe.get_nonzero_components(i)[comp_i] && 294 fe.get_nonzero_components(j)[comp_i]) 295 { 296 if (use_coefficient) 297 { 298 if (data.coefficient->n_components == 1) 299 for (unsigned int point = 0; point < n_q_points; 300 ++point) 301 add_data += 302 (data.coefficient_values[point] * 303 fe_values.shape_value_component(i, 304 point, 305 comp_i) * 306 fe_values.shape_value_component(j, 307 point, 308 comp_i) * 309 JxW[point]); 310 else 311 for (unsigned int point = 0; point < n_q_points; 312 ++point) 313 add_data += 314 (data.coefficient_vector_values[point](comp_i) * 315 fe_values.shape_value_component(i, 316 point, 317 comp_i) * 318 fe_values.shape_value_component(j, 319 point, 320 comp_i) * 321 JxW[point]); 322 } 323 else 324 for (unsigned int point = 0; point < n_q_points; 325 ++point) 326 add_data += 327 fe_values.shape_value_component(i, point, comp_i) * 328 fe_values.shape_value_component(j, point, comp_i) * 329 JxW[point]; 330 } 331 332 copy_data.cell_matrix(i, j) = add_data; 333 copy_data.cell_matrix(j, i) = add_data; 334 } 335 336 if (use_rhs_function) 337 { 338 number add_data = 0; 339 for (unsigned int comp_i = 0; comp_i < n_components; ++comp_i) 340 if (fe.get_nonzero_components(i)[comp_i]) 341 { 342 if (data.rhs_function->n_components == 1) 343 for (unsigned int point = 0; point < n_q_points; 344 ++point) 345 add_data += 346 data.rhs_values[point] * 347 fe_values.shape_value_component(i, point, comp_i) * 348 JxW[point]; 349 else 350 for (unsigned int point = 0; point < n_q_points; 351 ++point) 352 add_data += 353 data.rhs_vector_values[point](comp_i) * 354 fe_values.shape_value_component(i, point, comp_i) * 355 JxW[point]; 356 } 357 copy_data.cell_rhs(i) = add_data; 358 } 359 } 360 } 361 362 363 364 template <int dim, int spacedim, typename CellIterator> 365 void laplace_assembler(const CellIterator & cell,MatrixCreator::internal::AssemblerData::Scratch<dim,spacedim,double> & data,MatrixCreator::internal::AssemblerData::CopyData<double> & copy_data)366 laplace_assembler( 367 const CellIterator &cell, 368 MatrixCreator::internal::AssemblerData::Scratch<dim, spacedim, double> 369 & data, 370 MatrixCreator::internal::AssemblerData::CopyData<double> ©_data) 371 { 372 data.x_fe_values.reinit(cell); 373 const FEValues<dim, spacedim> &fe_values = 374 data.x_fe_values.get_present_fe_values(); 375 376 const unsigned int dofs_per_cell = fe_values.dofs_per_cell, 377 n_q_points = fe_values.n_quadrature_points; 378 const FiniteElement<dim, spacedim> &fe = fe_values.get_fe(); 379 const unsigned int n_components = fe.n_components(); 380 381 Assert(data.rhs_function == nullptr || 382 data.rhs_function->n_components == 1 || 383 data.rhs_function->n_components == n_components, 384 ::dealii::MatrixCreator::ExcComponentMismatch()); 385 Assert(data.coefficient == nullptr || 386 data.coefficient->n_components == 1 || 387 data.coefficient->n_components == n_components, 388 ::dealii::MatrixCreator::ExcComponentMismatch()); 389 390 copy_data.cell_matrix.reinit(dofs_per_cell, dofs_per_cell); 391 copy_data.cell_rhs.reinit(dofs_per_cell); 392 copy_data.dof_indices.resize(dofs_per_cell); 393 cell->get_dof_indices(copy_data.dof_indices); 394 395 396 const bool use_rhs_function = data.rhs_function != nullptr; 397 if (use_rhs_function) 398 { 399 if (data.rhs_function->n_components == 1) 400 { 401 data.rhs_values.resize(n_q_points); 402 data.rhs_function->value_list(fe_values.get_quadrature_points(), 403 data.rhs_values); 404 } 405 else 406 { 407 data.rhs_vector_values.resize( 408 n_q_points, dealii::Vector<double>(n_components)); 409 data.rhs_function->vector_value_list( 410 fe_values.get_quadrature_points(), data.rhs_vector_values); 411 } 412 } 413 414 const bool use_coefficient = data.coefficient != nullptr; 415 if (use_coefficient) 416 { 417 if (data.coefficient->n_components == 1) 418 { 419 data.coefficient_values.resize(n_q_points); 420 data.coefficient->value_list(fe_values.get_quadrature_points(), 421 data.coefficient_values); 422 } 423 else 424 { 425 data.coefficient_vector_values.resize( 426 n_q_points, dealii::Vector<double>(n_components)); 427 data.coefficient->vector_value_list( 428 fe_values.get_quadrature_points(), 429 data.coefficient_vector_values); 430 } 431 } 432 433 434 const std::vector<double> &JxW = fe_values.get_JxW_values(); 435 double add_data; 436 for (unsigned int i = 0; i < dofs_per_cell; ++i) 437 if (fe.is_primitive()) 438 { 439 const unsigned int component_i = 440 fe.system_to_component_index(i).first; 441 const Tensor<1, spacedim> *grad_phi_i = &fe_values.shape_grad(i, 0); 442 443 // can use symmetry 444 for (unsigned int j = i; j < dofs_per_cell; ++j) 445 if ((n_components == 1) || 446 (fe.system_to_component_index(j).first == component_i)) 447 { 448 const Tensor<1, spacedim> *grad_phi_j = 449 &fe_values.shape_grad(j, 0); 450 add_data = 0; 451 if (use_coefficient) 452 { 453 if (data.coefficient->n_components == 1) 454 for (unsigned int point = 0; point < n_q_points; 455 ++point) 456 add_data += 457 ((grad_phi_i[point] * grad_phi_j[point]) * 458 JxW[point] * data.coefficient_values[point]); 459 else 460 for (unsigned int point = 0; point < n_q_points; 461 ++point) 462 add_data += ((grad_phi_i[point] * grad_phi_j[point]) * 463 JxW[point] * 464 data.coefficient_vector_values[point]( 465 component_i)); 466 } 467 else 468 for (unsigned int point = 0; point < n_q_points; ++point) 469 add_data += 470 (grad_phi_i[point] * grad_phi_j[point]) * JxW[point]; 471 472 copy_data.cell_matrix(i, j) = add_data; 473 copy_data.cell_matrix(j, i) = add_data; 474 } 475 476 if (use_rhs_function) 477 { 478 const double *phi_i = &fe_values.shape_value(i, 0); 479 add_data = 0; 480 if (data.rhs_function->n_components == 1) 481 for (unsigned int point = 0; point < n_q_points; ++point) 482 add_data += 483 phi_i[point] * JxW[point] * data.rhs_values[point]; 484 else 485 for (unsigned int point = 0; point < n_q_points; ++point) 486 add_data += phi_i[point] * JxW[point] * 487 data.rhs_vector_values[point](component_i); 488 copy_data.cell_rhs(i) = add_data; 489 } 490 } 491 else 492 { 493 // non-primitive vector-valued FE 494 for (unsigned int j = i; j < dofs_per_cell; ++j) 495 { 496 add_data = 0; 497 for (unsigned int comp_i = 0; comp_i < n_components; ++comp_i) 498 if (fe.get_nonzero_components(i)[comp_i] && 499 fe.get_nonzero_components(j)[comp_i]) 500 { 501 if (use_coefficient) 502 { 503 if (data.coefficient->n_components == 1) 504 for (unsigned int point = 0; point < n_q_points; 505 ++point) 506 add_data += 507 ((fe_values.shape_grad_component(i, 508 point, 509 comp_i) * 510 fe_values.shape_grad_component(j, 511 point, 512 comp_i)) * 513 JxW[point] * data.coefficient_values[point]); 514 else 515 for (unsigned int point = 0; point < n_q_points; 516 ++point) 517 add_data += 518 ((fe_values.shape_grad_component(i, 519 point, 520 comp_i) * 521 fe_values.shape_grad_component(j, 522 point, 523 comp_i)) * 524 JxW[point] * 525 data.coefficient_vector_values[point](comp_i)); 526 } 527 else 528 for (unsigned int point = 0; point < n_q_points; 529 ++point) 530 add_data += 531 (fe_values.shape_grad_component(i, point, comp_i) * 532 fe_values.shape_grad_component(j, point, comp_i)) * 533 JxW[point]; 534 } 535 536 copy_data.cell_matrix(i, j) = add_data; 537 copy_data.cell_matrix(j, i) = add_data; 538 } 539 540 if (use_rhs_function) 541 { 542 add_data = 0; 543 for (unsigned int comp_i = 0; comp_i < n_components; ++comp_i) 544 if (fe.get_nonzero_components(i)[comp_i]) 545 { 546 if (data.rhs_function->n_components == 1) 547 for (unsigned int point = 0; point < n_q_points; 548 ++point) 549 add_data += 550 fe_values.shape_value_component(i, point, comp_i) * 551 JxW[point] * data.rhs_values[point]; 552 else 553 for (unsigned int point = 0; point < n_q_points; 554 ++point) 555 add_data += 556 fe_values.shape_value_component(i, point, comp_i) * 557 JxW[point] * data.rhs_vector_values[point](comp_i); 558 } 559 copy_data.cell_rhs(i) = add_data; 560 } 561 } 562 } 563 564 565 566 template <typename number, typename MatrixType, typename VectorType> 567 void copy_local_to_global(const AssemblerData::CopyData<number> & data,MatrixType * matrix,VectorType * right_hand_side)568 copy_local_to_global(const AssemblerData::CopyData<number> &data, 569 MatrixType * matrix, 570 VectorType * right_hand_side) 571 { 572 const unsigned int dofs_per_cell = data.dof_indices.size(); 573 (void)dofs_per_cell; 574 575 Assert(data.cell_matrix.m() == dofs_per_cell, ExcInternalError()); 576 Assert(data.cell_matrix.n() == dofs_per_cell, ExcInternalError()); 577 Assert((right_hand_side == nullptr) || 578 (data.cell_rhs.size() == dofs_per_cell), 579 ExcInternalError()); 580 581 if (right_hand_side != nullptr) 582 data.constraints->distribute_local_to_global(data.cell_matrix, 583 data.cell_rhs, 584 data.dof_indices, 585 *matrix, 586 *right_hand_side); 587 else 588 data.constraints->distribute_local_to_global(data.cell_matrix, 589 data.dof_indices, 590 *matrix); 591 } 592 593 594 595 namespace AssemblerBoundary 596 { 597 struct Scratch 598 { 599 Scratch() = default; 600 }; 601 602 template <int dim, int spacedim, typename number> 603 struct CopyData 604 { 605 CopyData(); 606 607 unsigned int dofs_per_cell; 608 std::vector<types::global_dof_index> dofs; 609 std::vector<std::vector<bool>> dof_is_on_face; 610 typename DoFHandler<dim, spacedim>::active_cell_iterator cell; 611 std::vector<FullMatrix<number>> cell_matrix; 612 std::vector<Vector<number>> cell_vector; 613 }; 614 615 616 template <int dim, int spacedim, typename number> CopyData()617 CopyData<dim, spacedim, number>::CopyData() 618 : dofs_per_cell(numbers::invalid_unsigned_int) 619 {} 620 621 622 } // namespace AssemblerBoundary 623 } // namespace internal 624 } // namespace MatrixCreator 625 626 627 namespace MatrixCreator 628 { 629 template <int dim, int spacedim, typename number> 630 void create_mass_matrix(const Mapping<dim,spacedim> & mapping,const DoFHandler<dim,spacedim> & dof,const Quadrature<dim> & q,SparseMatrix<number> & matrix,const Function<spacedim,number> * const coefficient,const AffineConstraints<number> & constraints)631 create_mass_matrix(const Mapping<dim, spacedim> & mapping, 632 const DoFHandler<dim, spacedim> & dof, 633 const Quadrature<dim> & q, 634 SparseMatrix<number> & matrix, 635 const Function<spacedim, number> *const coefficient, 636 const AffineConstraints<number> & constraints) 637 { 638 Assert(matrix.m() == dof.n_dofs(), 639 ExcDimensionMismatch(matrix.m(), dof.n_dofs())); 640 Assert(matrix.n() == dof.n_dofs(), 641 ExcDimensionMismatch(matrix.n(), dof.n_dofs())); 642 643 hp::FECollection<dim, spacedim> fe_collection(dof.get_fe()); 644 hp::QCollection<dim> q_collection(q); 645 hp::MappingCollection<dim, spacedim> mapping_collection(mapping); 646 MatrixCreator::internal::AssemblerData::Scratch<dim, spacedim, number> 647 assembler_data(fe_collection, 648 update_values | update_JxW_values | 649 (coefficient != nullptr ? update_quadrature_points : 650 UpdateFlags(0)), 651 coefficient, 652 /*rhs_function=*/nullptr, 653 q_collection, 654 mapping_collection); 655 656 MatrixCreator::internal::AssemblerData::CopyData<number> copy_data; 657 copy_data.cell_matrix.reinit( 658 assembler_data.fe_collection.max_dofs_per_cell(), 659 assembler_data.fe_collection.max_dofs_per_cell()); 660 copy_data.cell_rhs.reinit(assembler_data.fe_collection.max_dofs_per_cell()); 661 copy_data.dof_indices.resize( 662 assembler_data.fe_collection.max_dofs_per_cell()); 663 copy_data.constraints = &constraints; 664 665 WorkStream::run( 666 dof.begin_active(), 667 static_cast<typename DoFHandler<dim, spacedim>::active_cell_iterator>( 668 dof.end()), 669 &MatrixCreator::internal::mass_assembler< 670 dim, 671 spacedim, 672 typename DoFHandler<dim, spacedim>::active_cell_iterator, 673 number>, 674 [&matrix](const internal::AssemblerData::CopyData<number> &data) { 675 MatrixCreator::internal::copy_local_to_global( 676 data, &matrix, static_cast<Vector<number> *>(nullptr)); 677 }, 678 assembler_data, 679 copy_data); 680 } 681 682 683 684 template <int dim, int spacedim, typename number> 685 void create_mass_matrix(const DoFHandler<dim,spacedim> & dof,const Quadrature<dim> & q,SparseMatrix<number> & matrix,const Function<spacedim,number> * const coefficient,const AffineConstraints<number> & constraints)686 create_mass_matrix(const DoFHandler<dim, spacedim> & dof, 687 const Quadrature<dim> & q, 688 SparseMatrix<number> & matrix, 689 const Function<spacedim, number> *const coefficient, 690 const AffineConstraints<number> & constraints) 691 { 692 create_mass_matrix(StaticMappingQ1<dim, spacedim>::mapping, 693 dof, 694 q, 695 matrix, 696 coefficient, 697 constraints); 698 } 699 700 701 702 template <int dim, int spacedim, typename number> 703 void create_mass_matrix(const Mapping<dim,spacedim> & mapping,const DoFHandler<dim,spacedim> & dof,const Quadrature<dim> & q,SparseMatrix<number> & matrix,const Function<spacedim,number> & rhs,Vector<number> & rhs_vector,const Function<spacedim,number> * const coefficient,const AffineConstraints<number> & constraints)704 create_mass_matrix(const Mapping<dim, spacedim> & mapping, 705 const DoFHandler<dim, spacedim> & dof, 706 const Quadrature<dim> & q, 707 SparseMatrix<number> & matrix, 708 const Function<spacedim, number> & rhs, 709 Vector<number> & rhs_vector, 710 const Function<spacedim, number> *const coefficient, 711 const AffineConstraints<number> & constraints) 712 { 713 Assert(matrix.m() == dof.n_dofs(), 714 ExcDimensionMismatch(matrix.m(), dof.n_dofs())); 715 Assert(matrix.n() == dof.n_dofs(), 716 ExcDimensionMismatch(matrix.n(), dof.n_dofs())); 717 718 hp::FECollection<dim, spacedim> fe_collection(dof.get_fe()); 719 hp::QCollection<dim> q_collection(q); 720 hp::MappingCollection<dim, spacedim> mapping_collection(mapping); 721 MatrixCreator::internal::AssemblerData::Scratch<dim, spacedim, number> 722 assembler_data(fe_collection, 723 update_values | update_JxW_values | 724 update_quadrature_points, 725 coefficient, 726 &rhs, 727 q_collection, 728 mapping_collection); 729 MatrixCreator::internal::AssemblerData::CopyData<number> copy_data; 730 copy_data.cell_matrix.reinit( 731 assembler_data.fe_collection.max_dofs_per_cell(), 732 assembler_data.fe_collection.max_dofs_per_cell()); 733 copy_data.cell_rhs.reinit(assembler_data.fe_collection.max_dofs_per_cell()); 734 copy_data.dof_indices.resize( 735 assembler_data.fe_collection.max_dofs_per_cell()); 736 copy_data.constraints = &constraints; 737 738 WorkStream::run( 739 dof.begin_active(), 740 static_cast<typename DoFHandler<dim, spacedim>::active_cell_iterator>( 741 dof.end()), 742 &MatrixCreator::internal::mass_assembler< 743 dim, 744 spacedim, 745 typename DoFHandler<dim, spacedim>::active_cell_iterator, 746 number>, 747 [&matrix, 748 &rhs_vector](const internal::AssemblerData::CopyData<number> &data) { 749 MatrixCreator::internal::copy_local_to_global(data, 750 &matrix, 751 &rhs_vector); 752 }, 753 assembler_data, 754 copy_data); 755 } 756 757 758 759 template <int dim, int spacedim, typename number> 760 void create_mass_matrix(const DoFHandler<dim,spacedim> & dof,const Quadrature<dim> & q,SparseMatrix<number> & matrix,const Function<spacedim,number> & rhs,Vector<number> & rhs_vector,const Function<spacedim,number> * const coefficient,const AffineConstraints<number> & constraints)761 create_mass_matrix(const DoFHandler<dim, spacedim> & dof, 762 const Quadrature<dim> & q, 763 SparseMatrix<number> & matrix, 764 const Function<spacedim, number> & rhs, 765 Vector<number> & rhs_vector, 766 const Function<spacedim, number> *const coefficient, 767 const AffineConstraints<number> & constraints) 768 { 769 create_mass_matrix(StaticMappingQ1<dim, spacedim>::mapping, 770 dof, 771 q, 772 matrix, 773 rhs, 774 rhs_vector, 775 coefficient, 776 constraints); 777 } 778 779 780 781 template <int dim, int spacedim, typename number> 782 void create_mass_matrix(const hp::MappingCollection<dim,spacedim> & mapping,const DoFHandler<dim,spacedim> & dof,const hp::QCollection<dim> & q,SparseMatrix<number> & matrix,const Function<spacedim,number> * const coefficient,const AffineConstraints<number> & constraints)783 create_mass_matrix(const hp::MappingCollection<dim, spacedim> &mapping, 784 const DoFHandler<dim, spacedim> & dof, 785 const hp::QCollection<dim> & q, 786 SparseMatrix<number> & matrix, 787 const Function<spacedim, number> *const coefficient, 788 const AffineConstraints<number> & constraints) 789 { 790 Assert(matrix.m() == dof.n_dofs(), 791 ExcDimensionMismatch(matrix.m(), dof.n_dofs())); 792 Assert(matrix.n() == dof.n_dofs(), 793 ExcDimensionMismatch(matrix.n(), dof.n_dofs())); 794 795 MatrixCreator::internal::AssemblerData::Scratch<dim, spacedim, number> 796 assembler_data(dof.get_fe_collection(), 797 update_values | update_JxW_values | 798 (coefficient != nullptr ? update_quadrature_points : 799 UpdateFlags(0)), 800 coefficient, 801 /*rhs_function=*/nullptr, 802 q, 803 mapping); 804 MatrixCreator::internal::AssemblerData::CopyData<number> copy_data; 805 copy_data.cell_matrix.reinit( 806 assembler_data.fe_collection.max_dofs_per_cell(), 807 assembler_data.fe_collection.max_dofs_per_cell()); 808 copy_data.cell_rhs.reinit(assembler_data.fe_collection.max_dofs_per_cell()); 809 copy_data.dof_indices.resize( 810 assembler_data.fe_collection.max_dofs_per_cell()); 811 copy_data.constraints = &constraints; 812 813 WorkStream::run( 814 dof.begin_active(), 815 static_cast<typename DoFHandler<dim, spacedim>::active_cell_iterator>( 816 dof.end()), 817 &MatrixCreator::internal::mass_assembler< 818 dim, 819 spacedim, 820 typename DoFHandler<dim, spacedim>::active_cell_iterator, 821 number>, 822 [&matrix](const internal::AssemblerData::CopyData<number> &data) { 823 MatrixCreator::internal::copy_local_to_global( 824 data, &matrix, static_cast<Vector<number> *>(nullptr)); 825 }, 826 assembler_data, 827 copy_data); 828 } 829 830 831 832 template <int dim, int spacedim, typename number> 833 void create_mass_matrix(const DoFHandler<dim,spacedim> & dof,const hp::QCollection<dim> & q,SparseMatrix<number> & matrix,const Function<spacedim,number> * const coefficient,const AffineConstraints<number> & constraints)834 create_mass_matrix(const DoFHandler<dim, spacedim> & dof, 835 const hp::QCollection<dim> & q, 836 SparseMatrix<number> & matrix, 837 const Function<spacedim, number> *const coefficient, 838 const AffineConstraints<number> & constraints) 839 { 840 create_mass_matrix(hp::StaticMappingQ1<dim, spacedim>::mapping_collection, 841 dof, 842 q, 843 matrix, 844 coefficient, 845 constraints); 846 } 847 848 849 850 template <int dim, int spacedim, typename number> 851 void create_mass_matrix(const hp::MappingCollection<dim,spacedim> & mapping,const DoFHandler<dim,spacedim> & dof,const hp::QCollection<dim> & q,SparseMatrix<number> & matrix,const Function<spacedim,number> & rhs,Vector<number> & rhs_vector,const Function<spacedim,number> * const coefficient,const AffineConstraints<number> & constraints)852 create_mass_matrix(const hp::MappingCollection<dim, spacedim> &mapping, 853 const DoFHandler<dim, spacedim> & dof, 854 const hp::QCollection<dim> & q, 855 SparseMatrix<number> & matrix, 856 const Function<spacedim, number> & rhs, 857 Vector<number> & rhs_vector, 858 const Function<spacedim, number> *const coefficient, 859 const AffineConstraints<number> & constraints) 860 { 861 Assert(matrix.m() == dof.n_dofs(), 862 ExcDimensionMismatch(matrix.m(), dof.n_dofs())); 863 Assert(matrix.n() == dof.n_dofs(), 864 ExcDimensionMismatch(matrix.n(), dof.n_dofs())); 865 866 MatrixCreator::internal::AssemblerData::Scratch<dim, spacedim, number> 867 assembler_data(dof.get_fe_collection(), 868 update_values | update_JxW_values | 869 update_quadrature_points, 870 coefficient, 871 &rhs, 872 q, 873 mapping); 874 MatrixCreator::internal::AssemblerData::CopyData<number> copy_data; 875 copy_data.cell_matrix.reinit( 876 assembler_data.fe_collection.max_dofs_per_cell(), 877 assembler_data.fe_collection.max_dofs_per_cell()); 878 copy_data.cell_rhs.reinit(assembler_data.fe_collection.max_dofs_per_cell()); 879 copy_data.dof_indices.resize( 880 assembler_data.fe_collection.max_dofs_per_cell()); 881 copy_data.constraints = &constraints; 882 883 WorkStream::run( 884 dof.begin_active(), 885 static_cast<typename DoFHandler<dim, spacedim>::active_cell_iterator>( 886 dof.end()), 887 &MatrixCreator::internal::mass_assembler< 888 dim, 889 spacedim, 890 typename DoFHandler<dim, spacedim>::active_cell_iterator, 891 number>, 892 [&matrix, 893 &rhs_vector](const internal::AssemblerData::CopyData<number> &data) { 894 MatrixCreator::internal::copy_local_to_global(data, 895 &matrix, 896 &rhs_vector); 897 }, 898 assembler_data, 899 copy_data); 900 } 901 902 903 904 template <int dim, int spacedim, typename number> 905 void create_mass_matrix(const DoFHandler<dim,spacedim> & dof,const hp::QCollection<dim> & q,SparseMatrix<number> & matrix,const Function<spacedim,number> & rhs,Vector<number> & rhs_vector,const Function<spacedim,number> * const coefficient,const AffineConstraints<number> & constraints)906 create_mass_matrix(const DoFHandler<dim, spacedim> & dof, 907 const hp::QCollection<dim> & q, 908 SparseMatrix<number> & matrix, 909 const Function<spacedim, number> & rhs, 910 Vector<number> & rhs_vector, 911 const Function<spacedim, number> *const coefficient, 912 const AffineConstraints<number> & constraints) 913 { 914 create_mass_matrix(hp::StaticMappingQ1<dim, spacedim>::mapping_collection, 915 dof, 916 q, 917 matrix, 918 rhs, 919 rhs_vector, 920 coefficient, 921 constraints); 922 } 923 924 925 926 namespace internal 927 { 928 template <int dim, int spacedim, typename number> create_boundary_mass_matrix_1(typename DoFHandler<dim,spacedim>::active_cell_iterator const & cell,MatrixCreator::internal::AssemblerBoundary::Scratch const &,MatrixCreator::internal::AssemblerBoundary::CopyData<dim,spacedim,number> & copy_data,Mapping<dim,spacedim> const & mapping,FiniteElement<dim,spacedim> const & fe,Quadrature<dim-1> const & q,std::map<types::boundary_id,const Function<spacedim,number> * > const & boundary_functions,Function<spacedim,number> const * const coefficient,std::vector<unsigned int> const & component_mapping)929 void static inline create_boundary_mass_matrix_1( 930 typename DoFHandler<dim, spacedim>::active_cell_iterator const &cell, 931 MatrixCreator::internal::AssemblerBoundary::Scratch const &, 932 MatrixCreator::internal::AssemblerBoundary:: 933 CopyData<dim, spacedim, number> & copy_data, 934 Mapping<dim, spacedim> const & mapping, 935 FiniteElement<dim, spacedim> const &fe, 936 Quadrature<dim - 1> const & q, 937 std::map<types::boundary_id, const Function<spacedim, number> *> const 938 & boundary_functions, 939 Function<spacedim, number> const *const coefficient, 940 std::vector<unsigned int> const & component_mapping) 941 942 { 943 // Most assertions for this function are in the calling function 944 // before creating threads. 945 const unsigned int n_components = fe.n_components(); 946 const unsigned int n_function_components = 947 boundary_functions.begin()->second->n_components; 948 const bool fe_is_system = (n_components != 1); 949 const bool fe_is_primitive = fe.is_primitive(); 950 951 const unsigned int dofs_per_face = fe.n_dofs_per_face(); 952 953 copy_data.cell = cell; 954 copy_data.dofs_per_cell = fe.n_dofs_per_cell(); 955 956 UpdateFlags update_flags = 957 UpdateFlags(update_values | update_JxW_values | update_normal_vectors | 958 update_quadrature_points); 959 FEFaceValues<dim, spacedim> fe_values(mapping, fe, q, update_flags); 960 961 // two variables for the coefficient, one for the two cases 962 // indicated in the name 963 std::vector<number> coefficient_values(fe_values.n_quadrature_points, 1.); 964 std::vector<Vector<number>> coefficient_vector_values( 965 fe_values.n_quadrature_points, Vector<number>(n_components)); 966 const bool coefficient_is_vector = 967 (coefficient != nullptr && coefficient->n_components != 1); 968 969 std::vector<number> rhs_values_scalar(fe_values.n_quadrature_points); 970 std::vector<Vector<number>> rhs_values_system( 971 fe_values.n_quadrature_points, Vector<number>(n_function_components)); 972 973 copy_data.dofs.resize(copy_data.dofs_per_cell); 974 cell->get_dof_indices(copy_data.dofs); 975 976 std::vector<types::global_dof_index> dofs_on_face_vector(dofs_per_face); 977 978 // Because CopyData objects are reused and emplace_back is 979 // used, dof_is_on_face, cell_matrix, and cell_vector must be 980 // cleared before they are reused 981 copy_data.dof_is_on_face.clear(); 982 copy_data.cell_matrix.clear(); 983 copy_data.cell_vector.clear(); 984 985 for (const unsigned int face : cell->face_indices()) 986 // check if this face is on that part of the boundary we are 987 // interested in 988 if (boundary_functions.find(cell->face(face)->boundary_id()) != 989 boundary_functions.end()) 990 { 991 copy_data.cell_matrix.emplace_back(copy_data.dofs_per_cell, 992 copy_data.dofs_per_cell); 993 copy_data.cell_vector.emplace_back(copy_data.dofs_per_cell); 994 fe_values.reinit(cell, face); 995 996 if (fe_is_system) 997 // FE has several components 998 { 999 boundary_functions.find(cell->face(face)->boundary_id()) 1000 ->second->vector_value_list(fe_values.get_quadrature_points(), 1001 rhs_values_system); 1002 1003 if (coefficient_is_vector) 1004 // If coefficient is vector valued, fill all 1005 // components 1006 coefficient->vector_value_list( 1007 fe_values.get_quadrature_points(), 1008 coefficient_vector_values); 1009 else 1010 { 1011 // If a scalar function is given, update the 1012 // values, if not, use the default one set in the 1013 // constructor above 1014 if (coefficient != nullptr) 1015 coefficient->value_list(fe_values.get_quadrature_points(), 1016 coefficient_values); 1017 // Copy scalar values into vector 1018 for (unsigned int point = 0; 1019 point < fe_values.n_quadrature_points; 1020 ++point) 1021 coefficient_vector_values[point] = 1022 coefficient_values[point]; 1023 } 1024 1025 // Special treatment for Hdiv elements, 1026 // where only normal components should be projected. 1027 // For Hdiv we need to compute (u dot n, v dot n) which 1028 // can be done as sum over dim components of 1029 // u[c] * n[c] * v[c] * n[c] = u[c] * v[c] * 1030 // normal_adjustment[c] Same approach does not work for Hcurl, 1031 // so we throw an exception. Default value 1.0 allows for use 1032 // with non Hdiv elements 1033 std::vector<std::vector<double>> normal_adjustment( 1034 fe_values.n_quadrature_points, 1035 std::vector<double>(n_components, 1.)); 1036 1037 for (unsigned int comp = 0; comp < n_components; ++comp) 1038 { 1039 const FiniteElement<dim, spacedim> &base = 1040 fe.base_element(fe.component_to_base_index(comp).first); 1041 const unsigned int bcomp = 1042 fe.component_to_base_index(comp).second; 1043 1044 if (!base.conforms(FiniteElementData<dim>::H1) && 1045 base.conforms(FiniteElementData<dim>::Hdiv) && 1046 fe_is_primitive) 1047 Assert(false, ExcNotImplemented()); 1048 1049 if (!base.conforms(FiniteElementData<dim>::H1) && 1050 base.conforms(FiniteElementData<dim>::Hcurl)) 1051 Assert(false, ExcNotImplemented()); 1052 1053 if (!base.conforms(FiniteElementData<dim>::H1) && 1054 base.conforms(FiniteElementData<dim>::Hdiv)) 1055 for (unsigned int point = 0; 1056 point < fe_values.n_quadrature_points; 1057 ++point) 1058 normal_adjustment[point][comp] = 1059 fe_values.normal_vector(point)[bcomp] * 1060 fe_values.normal_vector(point)[bcomp]; 1061 } 1062 1063 for (unsigned int point = 0; 1064 point < fe_values.n_quadrature_points; 1065 ++point) 1066 { 1067 const double weight = fe_values.JxW(point); 1068 for (unsigned int i = 0; i < fe_values.dofs_per_cell; ++i) 1069 if (fe_is_primitive) 1070 { 1071 for (unsigned int j = 0; j < fe_values.dofs_per_cell; 1072 ++j) 1073 { 1074 if (fe.system_to_component_index(j).first == 1075 fe.system_to_component_index(i).first) 1076 { 1077 copy_data.cell_matrix.back()(i, j) += 1078 coefficient_vector_values[point]( 1079 fe.system_to_component_index(i).first) * 1080 weight * fe_values.shape_value(j, point) * 1081 fe_values.shape_value(i, point); 1082 } 1083 } 1084 copy_data.cell_vector.back()(i) += 1085 rhs_values_system[point]( 1086 component_mapping[fe.system_to_component_index(i) 1087 .first]) * 1088 fe_values.shape_value(i, point) * weight; 1089 } 1090 else 1091 { 1092 for (unsigned int comp = 0; comp < n_components; 1093 ++comp) 1094 { 1095 for (unsigned int j = 0; 1096 j < fe_values.dofs_per_cell; 1097 ++j) 1098 copy_data.cell_matrix.back()(i, j) += 1099 coefficient_vector_values[point](comp) * 1100 fe_values.shape_value_component(j, 1101 point, 1102 comp) * 1103 fe_values.shape_value_component(i, 1104 point, 1105 comp) * 1106 normal_adjustment[point][comp] * weight; 1107 copy_data.cell_vector.back()(i) += 1108 rhs_values_system[point]( 1109 component_mapping[comp]) * 1110 fe_values.shape_value_component(i, 1111 point, 1112 comp) * 1113 normal_adjustment[point][comp] * weight; 1114 } 1115 } 1116 } 1117 } 1118 else 1119 // FE is a scalar one 1120 { 1121 boundary_functions.find(cell->face(face)->boundary_id()) 1122 ->second->value_list(fe_values.get_quadrature_points(), 1123 rhs_values_scalar); 1124 1125 if (coefficient != nullptr) 1126 coefficient->value_list(fe_values.get_quadrature_points(), 1127 coefficient_values); 1128 for (unsigned int point = 0; 1129 point < fe_values.n_quadrature_points; 1130 ++point) 1131 { 1132 const double weight = fe_values.JxW(point); 1133 for (unsigned int i = 0; i < fe_values.dofs_per_cell; ++i) 1134 { 1135 const double v = fe_values.shape_value(i, point); 1136 for (unsigned int j = 0; j < fe_values.dofs_per_cell; 1137 ++j) 1138 { 1139 const double u = fe_values.shape_value(j, point); 1140 copy_data.cell_matrix.back()(i, j) += 1141 (coefficient_values[point] * u * v * weight); 1142 } 1143 copy_data.cell_vector.back()(i) += 1144 rhs_values_scalar[point] * v * weight; 1145 } 1146 } 1147 } 1148 1149 1150 cell->face(face)->get_dof_indices(dofs_on_face_vector); 1151 // for each dof on the cell, have a flag whether it is on 1152 // the face 1153 copy_data.dof_is_on_face.emplace_back(copy_data.dofs_per_cell); 1154 // check for each of the dofs on this cell whether it is 1155 // on the face 1156 for (unsigned int i = 0; i < copy_data.dofs_per_cell; ++i) 1157 copy_data.dof_is_on_face.back()[i] = 1158 (std::find(dofs_on_face_vector.begin(), 1159 dofs_on_face_vector.end(), 1160 copy_data.dofs[i]) != dofs_on_face_vector.end()); 1161 } 1162 } 1163 1164 template <int dim, int spacedim, typename number> 1165 void copy_boundary_mass_matrix_1(MatrixCreator::internal::AssemblerBoundary::CopyData<dim,spacedim,number> const & copy_data,std::map<types::boundary_id,const Function<spacedim,number> * > const & boundary_functions,std::vector<types::global_dof_index> const & dof_to_boundary_mapping,SparseMatrix<number> & matrix,Vector<number> & rhs_vector)1166 copy_boundary_mass_matrix_1( 1167 MatrixCreator::internal::AssemblerBoundary:: 1168 CopyData<dim, spacedim, number> const ©_data, 1169 std::map<types::boundary_id, const Function<spacedim, number> *> const 1170 & boundary_functions, 1171 std::vector<types::global_dof_index> const &dof_to_boundary_mapping, 1172 SparseMatrix<number> & matrix, 1173 Vector<number> & rhs_vector) 1174 { 1175 // now transfer cell matrix and vector to the whole boundary matrix 1176 // 1177 // in the following: dof[i] holds the global index of the i-th degree of 1178 // freedom on the present cell. If it is also a dof on the boundary, it 1179 // must be a nonzero entry in the dof_to_boundary_mapping and then 1180 // the boundary index of this dof is dof_to_boundary_mapping[dof[i]]. 1181 // 1182 // if dof[i] is not on the boundary, it should be zero on the boundary 1183 // therefore on all quadrature points and finally all of its 1184 // entries in the cell matrix and vector should be zero. If not, we 1185 // throw an error (note: because of the evaluation of the shape 1186 // functions only up to machine precision, the term "must be zero" 1187 // really should mean: "should be very small". since this is only an 1188 // assertion and not part of the code, we may choose "very small" 1189 // quite arbitrarily) 1190 // 1191 // the main problem here is that the matrix or vector entry should also 1192 // be zero if the degree of freedom dof[i] is on the boundary, but not 1193 // on the present face, i.e. on another face of the same cell also 1194 // on the boundary. We can therefore not rely on the 1195 // dof_to_boundary_mapping[dof[i]] being !=-1, we really have to 1196 // determine whether dof[i] is a dof on the present face. We do so 1197 // by getting the dofs on the face into @p{dofs_on_face_vector}, 1198 // a vector as always. Usually, searching in a vector is 1199 // inefficient, so we copy the dofs into a set, which enables binary 1200 // searches. 1201 unsigned int pos(0); 1202 for (const unsigned int face : copy_data.cell->face_indices()) 1203 { 1204 // check if this face is on that part of 1205 // the boundary we are interested in 1206 if (boundary_functions.find( 1207 copy_data.cell->face(face)->boundary_id()) != 1208 boundary_functions.end()) 1209 { 1210 for (unsigned int i = 0; i < copy_data.dofs_per_cell; ++i) 1211 { 1212 if (copy_data.dof_is_on_face[pos][i] && 1213 dof_to_boundary_mapping[copy_data.dofs[i]] != 1214 numbers::invalid_dof_index) 1215 { 1216 for (unsigned int j = 0; j < copy_data.dofs_per_cell; ++j) 1217 if (copy_data.dof_is_on_face[pos][j] && 1218 dof_to_boundary_mapping[copy_data.dofs[j]] != 1219 numbers::invalid_dof_index) 1220 { 1221 AssertIsFinite(copy_data.cell_matrix[pos](i, j)); 1222 matrix.add( 1223 dof_to_boundary_mapping[copy_data.dofs[i]], 1224 dof_to_boundary_mapping[copy_data.dofs[j]], 1225 copy_data.cell_matrix[pos](i, j)); 1226 } 1227 AssertIsFinite(copy_data.cell_vector[pos](i)); 1228 rhs_vector(dof_to_boundary_mapping[copy_data.dofs[i]]) += 1229 copy_data.cell_vector[pos](i); 1230 } 1231 } 1232 ++pos; 1233 } 1234 } 1235 } 1236 1237 1238 template <> 1239 void inline create_boundary_mass_matrix_1<1, 3, float>( 1240 DoFHandler<1, 3>::active_cell_iterator const & /*cell*/, 1241 MatrixCreator::internal::AssemblerBoundary::Scratch const &, 1242 MatrixCreator::internal::AssemblerBoundary::CopyData<1, 3, float> 1243 & /*copy_data*/, 1244 Mapping<1, 3> const &, 1245 FiniteElement<1, 3> const &, 1246 Quadrature<0> const &, 1247 std::map<types::boundary_id, const Function<3, float> *> const 1248 & /*boundary_functions*/, 1249 Function<3, float> const *const /*coefficient*/, 1250 std::vector<unsigned int> const & /*component_mapping*/) 1251 { 1252 Assert(false, ExcNotImplemented()); 1253 } 1254 1255 template <> 1256 void inline create_boundary_mass_matrix_1<1, 3, double>( 1257 DoFHandler<1, 3>::active_cell_iterator const & /*cell*/, 1258 MatrixCreator::internal::AssemblerBoundary::Scratch const &, 1259 MatrixCreator::internal::AssemblerBoundary::CopyData<1, 3, double> 1260 & /*copy_data*/, 1261 Mapping<1, 3> const &, 1262 FiniteElement<1, 3> const &, 1263 Quadrature<0> const &, 1264 std::map<types::boundary_id, const Function<3, double> *> const 1265 & /*boundary_functions*/, 1266 Function<3, double> const *const /*coefficient*/, 1267 std::vector<unsigned int> const & /*component_mapping*/) 1268 { 1269 Assert(false, ExcNotImplemented()); 1270 } 1271 1272 } // namespace internal 1273 1274 1275 1276 template <int dim, int spacedim, typename number> 1277 void create_boundary_mass_matrix(const Mapping<dim,spacedim> & mapping,const DoFHandler<dim,spacedim> & dof,const Quadrature<dim-1> & q,SparseMatrix<number> & matrix,const std::map<types::boundary_id,const Function<spacedim,number> * > & boundary_functions,Vector<number> & rhs_vector,std::vector<types::global_dof_index> & dof_to_boundary_mapping,const Function<spacedim,number> * const coefficient,std::vector<unsigned int> component_mapping)1278 create_boundary_mass_matrix( 1279 const Mapping<dim, spacedim> & mapping, 1280 const DoFHandler<dim, spacedim> &dof, 1281 const Quadrature<dim - 1> & q, 1282 SparseMatrix<number> & matrix, 1283 const std::map<types::boundary_id, const Function<spacedim, number> *> 1284 & boundary_functions, 1285 Vector<number> & rhs_vector, 1286 std::vector<types::global_dof_index> & dof_to_boundary_mapping, 1287 const Function<spacedim, number> *const coefficient, 1288 std::vector<unsigned int> component_mapping) 1289 { 1290 // what would that be in 1d? the 1291 // identity matrix on the boundary 1292 // dofs? 1293 if (dim == 1) 1294 { 1295 Assert(false, ExcNotImplemented()); 1296 return; 1297 } 1298 1299 const FiniteElement<dim, spacedim> &fe = dof.get_fe(); 1300 const unsigned int n_components = fe.n_components(); 1301 1302 Assert(matrix.n() == dof.n_boundary_dofs(boundary_functions), 1303 ExcInternalError()); 1304 Assert(matrix.n() == matrix.m(), ExcInternalError()); 1305 Assert(matrix.n() == rhs_vector.size(), ExcInternalError()); 1306 Assert(boundary_functions.size() != 0, ExcInternalError()); 1307 Assert(dof_to_boundary_mapping.size() == dof.n_dofs(), ExcInternalError()); 1308 Assert(coefficient == nullptr || coefficient->n_components == 1 || 1309 coefficient->n_components == n_components, 1310 ExcComponentMismatch()); 1311 1312 if (component_mapping.size() == 0) 1313 { 1314 AssertDimension(n_components, 1315 boundary_functions.begin()->second->n_components); 1316 for (unsigned int i = 0; i < n_components; ++i) 1317 component_mapping.push_back(i); 1318 } 1319 else 1320 AssertDimension(n_components, component_mapping.size()); 1321 1322 MatrixCreator::internal::AssemblerBoundary::Scratch scratch; 1323 MatrixCreator::internal::AssemblerBoundary::CopyData<dim, spacedim, number> 1324 copy_data; 1325 1326 WorkStream::run( 1327 dof.begin_active(), 1328 dof.end(), 1329 [&mapping, &fe, &q, &boundary_functions, coefficient, &component_mapping]( 1330 typename DoFHandler<dim, spacedim>::active_cell_iterator const &cell, 1331 MatrixCreator::internal::AssemblerBoundary::Scratch const &scratch_data, 1332 MatrixCreator::internal::AssemblerBoundary:: 1333 CopyData<dim, spacedim, number> ©_data) { 1334 internal::create_boundary_mass_matrix_1(cell, 1335 scratch_data, 1336 copy_data, 1337 mapping, 1338 fe, 1339 q, 1340 boundary_functions, 1341 coefficient, 1342 component_mapping); 1343 }, 1344 [&boundary_functions, &dof_to_boundary_mapping, &matrix, &rhs_vector]( 1345 MatrixCreator::internal::AssemblerBoundary:: 1346 CopyData<dim, spacedim, number> const ©_data) { 1347 internal::copy_boundary_mass_matrix_1(copy_data, 1348 boundary_functions, 1349 dof_to_boundary_mapping, 1350 matrix, 1351 rhs_vector); 1352 }, 1353 scratch, 1354 copy_data); 1355 } 1356 1357 1358 1359 namespace internal 1360 { 1361 template <int dim, int spacedim, typename number> 1362 void create_hp_boundary_mass_matrix_1(typename DoFHandler<dim,spacedim>::active_cell_iterator const & cell,MatrixCreator::internal::AssemblerBoundary::Scratch const &,MatrixCreator::internal::AssemblerBoundary::CopyData<dim,spacedim,number> & copy_data,hp::MappingCollection<dim,spacedim> const & mapping,hp::FECollection<dim,spacedim> const & fe_collection,hp::QCollection<dim-1> const & q,const std::map<types::boundary_id,const Function<spacedim,number> * > & boundary_functions,Function<spacedim,number> const * const coefficient,std::vector<unsigned int> const & component_mapping)1363 create_hp_boundary_mass_matrix_1( 1364 typename DoFHandler<dim, spacedim>::active_cell_iterator const &cell, 1365 MatrixCreator::internal::AssemblerBoundary::Scratch const &, 1366 MatrixCreator::internal::AssemblerBoundary :: 1367 CopyData<dim, spacedim, number> & copy_data, 1368 hp::MappingCollection<dim, spacedim> const &mapping, 1369 hp::FECollection<dim, spacedim> const & fe_collection, 1370 hp::QCollection<dim - 1> const & q, 1371 const std::map<types::boundary_id, const Function<spacedim, number> *> 1372 & boundary_functions, 1373 Function<spacedim, number> const *const coefficient, 1374 std::vector<unsigned int> const & component_mapping) 1375 { 1376 const unsigned int n_components = fe_collection.n_components(); 1377 const unsigned int n_function_components = 1378 boundary_functions.begin()->second->n_components; 1379 const FiniteElement<dim, spacedim> &fe = cell->get_fe(); 1380 const bool fe_is_system = (n_components != 1); 1381 const bool fe_is_primitive = fe.is_primitive(); 1382 const unsigned int dofs_per_face = fe.n_dofs_per_face(); 1383 1384 copy_data.cell = cell; 1385 copy_data.dofs_per_cell = fe.n_dofs_per_cell(); 1386 copy_data.dofs.resize(copy_data.dofs_per_cell); 1387 cell->get_dof_indices(copy_data.dofs); 1388 1389 1390 UpdateFlags update_flags = 1391 UpdateFlags(update_values | update_JxW_values | update_normal_vectors | 1392 update_quadrature_points); 1393 hp::FEFaceValues<dim, spacedim> x_fe_values(mapping, 1394 fe_collection, 1395 q, 1396 update_flags); 1397 1398 // two variables for the coefficient, 1399 // one for the two cases indicated in 1400 // the name 1401 std::vector<number> coefficient_values; 1402 std::vector<Vector<number>> coefficient_vector_values; 1403 1404 const bool coefficient_is_vector = 1405 (coefficient != nullptr && coefficient->n_components != 1); 1406 1407 std::vector<number> rhs_values_scalar; 1408 std::vector<Vector<number>> rhs_values_system; 1409 1410 std::vector<types::global_dof_index> dofs_on_face_vector(dofs_per_face); 1411 1412 copy_data.dofs.resize(copy_data.dofs_per_cell); 1413 cell->get_dof_indices(copy_data.dofs); 1414 1415 // Because CopyData objects are reused and that push_back is used, 1416 // dof_is_on_face, cell_matrix, and cell_vector must be cleared before 1417 // they are reused 1418 copy_data.dof_is_on_face.clear(); 1419 copy_data.cell_matrix.clear(); 1420 copy_data.cell_vector.clear(); 1421 1422 1423 for (const unsigned int face : cell->face_indices()) 1424 // check if this face is on that part of 1425 // the boundary we are interested in 1426 if (boundary_functions.find(cell->face(face)->boundary_id()) != 1427 boundary_functions.end()) 1428 { 1429 x_fe_values.reinit(cell, face); 1430 1431 const FEFaceValues<dim, spacedim> &fe_values = 1432 x_fe_values.get_present_fe_values(); 1433 1434 copy_data.cell_matrix.emplace_back(copy_data.dofs_per_cell, 1435 copy_data.dofs_per_cell); 1436 copy_data.cell_vector.emplace_back(copy_data.dofs_per_cell); 1437 1438 if (fe_is_system) 1439 // FE has several components 1440 { 1441 coefficient_vector_values.resize(fe_values.n_quadrature_points, 1442 Vector<number>(n_components)); 1443 coefficient_values.resize(fe_values.n_quadrature_points, 1.); 1444 1445 rhs_values_system.resize(fe_values.n_quadrature_points, 1446 Vector<number>(n_function_components)); 1447 boundary_functions.find(cell->face(face)->boundary_id()) 1448 ->second->vector_value_list(fe_values.get_quadrature_points(), 1449 rhs_values_system); 1450 if (coefficient_is_vector) 1451 // In case coefficient is vector-valued, fill 1452 // all components 1453 coefficient->vector_value_list( 1454 fe_values.get_quadrature_points(), 1455 coefficient_vector_values); 1456 else 1457 // In case the scalar function is given, update the 1458 // values, if not - use the default (1.0) 1459 { 1460 if (coefficient != nullptr) 1461 coefficient->value_list(fe_values.get_quadrature_points(), 1462 coefficient_values); 1463 1464 for (unsigned int point = 0; 1465 point < fe_values.n_quadrature_points; 1466 ++point) 1467 coefficient_vector_values[point] = 1468 coefficient_values[point]; 1469 } 1470 1471 // Special treatment for Hdiv elements, 1472 // where only normal components should be projected. 1473 // For Hdiv we need to compute (u dot n, v dot n) which 1474 // can be done as sum over dim components of 1475 // u[c] * n[c] * v[c] * n[c] = u[c] * v[c] * 1476 // normal_adjustment[c] Same approach does not work for Hcurl, 1477 // so we throw an exception. Default value 1.0 allows for use 1478 // with non Hdiv elements 1479 std::vector<std::vector<double>> normal_adjustment( 1480 fe_values.n_quadrature_points, 1481 std::vector<double>(n_components, 1.)); 1482 1483 for (unsigned int comp = 0; comp < n_components; ++comp) 1484 { 1485 const FiniteElement<dim, spacedim> &base = 1486 fe.base_element(fe.component_to_base_index(comp).first); 1487 const unsigned int bcomp = 1488 fe.component_to_base_index(comp).second; 1489 1490 if (!base.conforms(FiniteElementData<dim>::H1) && 1491 base.conforms(FiniteElementData<dim>::Hdiv) && 1492 fe_is_primitive) 1493 Assert(false, ExcNotImplemented()); 1494 1495 if (!base.conforms(FiniteElementData<dim>::H1) && 1496 base.conforms(FiniteElementData<dim>::Hcurl)) 1497 Assert(false, ExcNotImplemented()); 1498 1499 if (!base.conforms(FiniteElementData<dim>::H1) && 1500 base.conforms(FiniteElementData<dim>::Hdiv)) 1501 for (unsigned int point = 0; 1502 point < fe_values.n_quadrature_points; 1503 ++point) 1504 normal_adjustment[point][comp] = 1505 fe_values.normal_vector(point)[bcomp] * 1506 fe_values.normal_vector(point)[bcomp]; 1507 } 1508 1509 for (unsigned int point = 0; 1510 point < fe_values.n_quadrature_points; 1511 ++point) 1512 { 1513 const double weight = fe_values.JxW(point); 1514 for (unsigned int i = 0; i < fe_values.dofs_per_cell; ++i) 1515 if (fe_is_primitive) 1516 { 1517 for (unsigned int j = 0; j < fe_values.dofs_per_cell; 1518 ++j) 1519 { 1520 if (fe.system_to_component_index(i).first == 1521 fe.system_to_component_index(j).first) 1522 { 1523 copy_data.cell_matrix.back()(i, j) += 1524 coefficient_vector_values[point]( 1525 fe.system_to_component_index(i).first) * 1526 fe_values.shape_value(i, point) * 1527 fe_values.shape_value(j, point) * weight; 1528 } 1529 } 1530 1531 copy_data.cell_vector.back()(i) += 1532 rhs_values_system[point]( 1533 component_mapping[fe.system_to_component_index(i) 1534 .first]) * 1535 fe_values.shape_value(i, point) * weight; 1536 } 1537 else 1538 { 1539 for (unsigned int comp = 0; comp < n_components; 1540 ++comp) 1541 { 1542 for (unsigned int j = 0; 1543 j < fe_values.dofs_per_cell; 1544 ++j) 1545 copy_data.cell_matrix.back()(i, j) += 1546 coefficient_vector_values[point](comp) * 1547 fe_values.shape_value_component(i, 1548 point, 1549 comp) * 1550 fe_values.shape_value_component(j, 1551 point, 1552 comp) * 1553 normal_adjustment[point][comp] * weight; 1554 copy_data.cell_vector.back()(i) += 1555 rhs_values_system[point]( 1556 component_mapping[comp]) * 1557 fe_values.shape_value_component(i, 1558 point, 1559 comp) * 1560 normal_adjustment[point][comp] * weight; 1561 } 1562 } 1563 } 1564 } 1565 else 1566 // FE is scalar 1567 { 1568 coefficient_values.resize(fe_values.n_quadrature_points, 1.); 1569 rhs_values_scalar.resize(fe_values.n_quadrature_points); 1570 boundary_functions.find(cell->face(face)->boundary_id()) 1571 ->second->value_list(fe_values.get_quadrature_points(), 1572 rhs_values_scalar); 1573 1574 if (coefficient != nullptr) 1575 { 1576 coefficient_values.resize(fe_values.n_quadrature_points); 1577 coefficient->value_list(fe_values.get_quadrature_points(), 1578 coefficient_values); 1579 } 1580 1581 for (unsigned int point = 0; 1582 point < fe_values.n_quadrature_points; 1583 ++point) 1584 { 1585 const double weight = fe_values.JxW(point); 1586 for (unsigned int i = 0; i < fe_values.dofs_per_cell; ++i) 1587 { 1588 const double v = fe_values.shape_value(i, point); 1589 for (unsigned int j = 0; j < fe_values.dofs_per_cell; 1590 ++j) 1591 { 1592 const double u = fe_values.shape_value(j, point); 1593 copy_data.cell_matrix.back()(i, j) += 1594 (coefficient_values[point] * v * u * weight); 1595 } 1596 copy_data.cell_vector.back()(i) += 1597 rhs_values_scalar[point] * v * weight; 1598 } 1599 } 1600 } 1601 1602 cell->face(face)->get_dof_indices(dofs_on_face_vector, 1603 cell->active_fe_index()); 1604 // for each dof on the cell, have a 1605 // flag whether it is on the face 1606 copy_data.dof_is_on_face.emplace_back(copy_data.dofs_per_cell); 1607 // check for each of the dofs on this cell 1608 // whether it is on the face 1609 for (unsigned int i = 0; i < copy_data.dofs_per_cell; ++i) 1610 copy_data.dof_is_on_face.back()[i] = 1611 (std::find(dofs_on_face_vector.begin(), 1612 dofs_on_face_vector.end(), 1613 copy_data.dofs[i]) != dofs_on_face_vector.end()); 1614 } 1615 } 1616 1617 1618 template <int dim, int spacedim, typename number> 1619 void copy_hp_boundary_mass_matrix_1(MatrixCreator::internal::AssemblerBoundary::CopyData<dim,spacedim,number> const & copy_data,std::map<types::boundary_id,const Function<spacedim,number> * > const & boundary_functions,std::vector<types::global_dof_index> const & dof_to_boundary_mapping,SparseMatrix<number> & matrix,Vector<number> & rhs_vector)1620 copy_hp_boundary_mass_matrix_1( 1621 MatrixCreator::internal::AssemblerBoundary :: 1622 CopyData<dim, spacedim, number> const ©_data, 1623 std::map<types::boundary_id, const Function<spacedim, number> *> const 1624 & boundary_functions, 1625 std::vector<types::global_dof_index> const &dof_to_boundary_mapping, 1626 SparseMatrix<number> & matrix, 1627 Vector<number> & rhs_vector) 1628 { 1629 // now transfer cell matrix and vector to the whole boundary matrix 1630 // 1631 // in the following: dof[i] holds the global index of the i-th degree of 1632 // freedom on the present cell. If it is also a dof on the boundary, it 1633 // must be a nonzero entry in the dof_to_boundary_mapping and then 1634 // the boundary index of this dof is dof_to_boundary_mapping[dof[i]]. 1635 // 1636 // if dof[i] is not on the boundary, it should be zero on the boundary 1637 // therefore on all quadrature points and finally all of its 1638 // entries in the cell matrix and vector should be zero. If not, we 1639 // throw an error (note: because of the evaluation of the shape 1640 // functions only up to machine precision, the term "must be zero" 1641 // really should mean: "should be very small". since this is only an 1642 // assertion and not part of the code, we may choose "very small" 1643 // quite arbitrarily) 1644 // 1645 // the main problem here is that the matrix or vector entry should also 1646 // be zero if the degree of freedom dof[i] is on the boundary, but not 1647 // on the present face, i.e. on another face of the same cell also 1648 // on the boundary. We can therefore not rely on the 1649 // dof_to_boundary_mapping[dof[i]] being !=-1, we really have to 1650 // determine whether dof[i] is a dof on the present face. We do so 1651 // by getting the dofs on the face into @p{dofs_on_face_vector}, 1652 // a vector as always. Usually, searching in a vector is 1653 // inefficient, so we copy the dofs into a set, which enables binary 1654 // searches. 1655 unsigned int pos(0); 1656 for (const unsigned int face : copy_data.cell->face_indices()) 1657 { 1658 // check if this face is on that part of 1659 // the boundary we are interested in 1660 if (boundary_functions.find( 1661 copy_data.cell->face(face)->boundary_id()) != 1662 boundary_functions.end()) 1663 { 1664 #ifdef DEBUG 1665 // in debug mode: compute an element in the matrix which is 1666 // guaranteed to belong to a boundary dof. We do this to check 1667 // that the entries in the cell matrix are guaranteed to be zero 1668 // if the respective dof is not on the boundary. Since because of 1669 // round-off, the actual value of the matrix entry may be 1670 // only close to zero, we assert that it is small relative to an 1671 // element which is guaranteed to be nonzero. (absolute smallness 1672 // does not suffice since the size of the domain scales in here) 1673 // 1674 // for this purpose we seek the diagonal of the matrix, where 1675 // there must be an element belonging to the boundary. we take the 1676 // maximum diagonal entry. 1677 types::global_dof_index max_element = 0; 1678 for (const auto index : dof_to_boundary_mapping) 1679 if ((index != numbers::invalid_dof_index) && 1680 (index > max_element)) 1681 max_element = index; 1682 Assert(max_element == matrix.n() - 1, ExcInternalError()); 1683 1684 double max_diag_entry = 0; 1685 for (unsigned int i = 0; i < copy_data.dofs_per_cell; ++i) 1686 if (std::abs(copy_data.cell_matrix[pos](i, i)) > max_diag_entry) 1687 max_diag_entry = std::abs(copy_data.cell_matrix[pos](i, i)); 1688 #endif 1689 1690 for (unsigned int i = 0; i < copy_data.dofs_per_cell; ++i) 1691 { 1692 if (copy_data.dof_is_on_face[pos][i] && 1693 dof_to_boundary_mapping[copy_data.dofs[i]] != 1694 numbers::invalid_dof_index) 1695 { 1696 for (unsigned int j = 0; j < copy_data.dofs_per_cell; ++j) 1697 if (copy_data.dof_is_on_face[pos][j] && 1698 dof_to_boundary_mapping[copy_data.dofs[j]] != 1699 numbers::invalid_dof_index) 1700 { 1701 AssertIsFinite(copy_data.cell_matrix[pos](i, j)); 1702 matrix.add( 1703 dof_to_boundary_mapping[copy_data.dofs[i]], 1704 dof_to_boundary_mapping[copy_data.dofs[j]], 1705 copy_data.cell_matrix[pos](i, j)); 1706 } 1707 AssertIsFinite(copy_data.cell_vector[pos](i)); 1708 rhs_vector(dof_to_boundary_mapping[copy_data.dofs[i]]) += 1709 copy_data.cell_vector[pos](i); 1710 } 1711 } 1712 ++pos; 1713 } 1714 } 1715 } 1716 } // namespace internal 1717 1718 1719 1720 template <int dim, int spacedim, typename number> 1721 void create_boundary_mass_matrix(const DoFHandler<dim,spacedim> & dof,const Quadrature<dim-1> & q,SparseMatrix<number> & matrix,const std::map<types::boundary_id,const Function<spacedim,number> * > & rhs,Vector<number> & rhs_vector,std::vector<types::global_dof_index> & dof_to_boundary_mapping,const Function<spacedim,number> * const a,std::vector<unsigned int> component_mapping)1722 create_boundary_mass_matrix( 1723 const DoFHandler<dim, spacedim> &dof, 1724 const Quadrature<dim - 1> & q, 1725 SparseMatrix<number> & matrix, 1726 const std::map<types::boundary_id, const Function<spacedim, number> *> &rhs, 1727 Vector<number> & rhs_vector, 1728 std::vector<types::global_dof_index> & dof_to_boundary_mapping, 1729 const Function<spacedim, number> *const a, 1730 std::vector<unsigned int> component_mapping) 1731 { 1732 create_boundary_mass_matrix(StaticMappingQ1<dim, spacedim>::mapping, 1733 dof, 1734 q, 1735 matrix, 1736 rhs, 1737 rhs_vector, 1738 dof_to_boundary_mapping, 1739 a, 1740 component_mapping); 1741 } 1742 1743 1744 1745 template <int dim, int spacedim, typename number> 1746 void create_boundary_mass_matrix(const hp::MappingCollection<dim,spacedim> & mapping,const DoFHandler<dim,spacedim> & dof,const hp::QCollection<dim-1> & q,SparseMatrix<number> & matrix,const std::map<types::boundary_id,const Function<spacedim,number> * > & boundary_functions,Vector<number> & rhs_vector,std::vector<types::global_dof_index> & dof_to_boundary_mapping,const Function<spacedim,number> * const coefficient,std::vector<unsigned int> component_mapping)1747 create_boundary_mass_matrix( 1748 const hp::MappingCollection<dim, spacedim> &mapping, 1749 const DoFHandler<dim, spacedim> & dof, 1750 const hp::QCollection<dim - 1> & q, 1751 SparseMatrix<number> & matrix, 1752 const std::map<types::boundary_id, const Function<spacedim, number> *> 1753 & boundary_functions, 1754 Vector<number> & rhs_vector, 1755 std::vector<types::global_dof_index> & dof_to_boundary_mapping, 1756 const Function<spacedim, number> *const coefficient, 1757 std::vector<unsigned int> component_mapping) 1758 { 1759 // what would that be in 1d? the 1760 // identity matrix on the boundary 1761 // dofs? 1762 if (dim == 1) 1763 { 1764 Assert(false, ExcNotImplemented()); 1765 return; 1766 } 1767 1768 const hp::FECollection<dim, spacedim> &fe_collection = 1769 dof.get_fe_collection(); 1770 const unsigned int n_components = fe_collection.n_components(); 1771 1772 Assert(matrix.n() == dof.n_boundary_dofs(boundary_functions), 1773 ExcInternalError()); 1774 Assert(matrix.n() == matrix.m(), ExcInternalError()); 1775 Assert(matrix.n() == rhs_vector.size(), ExcInternalError()); 1776 Assert(boundary_functions.size() != 0, ExcInternalError()); 1777 Assert(dof_to_boundary_mapping.size() == dof.n_dofs(), ExcInternalError()); 1778 Assert(coefficient == nullptr || coefficient->n_components == 1 || 1779 coefficient->n_components == n_components, 1780 ExcComponentMismatch()); 1781 1782 if (component_mapping.size() == 0) 1783 { 1784 AssertDimension(n_components, 1785 boundary_functions.begin()->second->n_components); 1786 for (unsigned int i = 0; i < n_components; ++i) 1787 component_mapping.push_back(i); 1788 } 1789 else 1790 AssertDimension(n_components, component_mapping.size()); 1791 1792 MatrixCreator::internal::AssemblerBoundary::Scratch scratch; 1793 MatrixCreator::internal::AssemblerBoundary::CopyData<dim, spacedim, number> 1794 copy_data; 1795 1796 WorkStream::run( 1797 dof.begin_active(), 1798 dof.end(), 1799 [&mapping, 1800 &fe_collection, 1801 &q, 1802 &boundary_functions, 1803 coefficient, 1804 &component_mapping]( 1805 typename DoFHandler<dim, spacedim>::active_cell_iterator const &cell, 1806 MatrixCreator::internal::AssemblerBoundary::Scratch const &scratch_data, 1807 MatrixCreator::internal::AssemblerBoundary :: 1808 CopyData<dim, spacedim, number> ©_data) { 1809 internal::create_hp_boundary_mass_matrix_1(cell, 1810 scratch_data, 1811 copy_data, 1812 mapping, 1813 fe_collection, 1814 q, 1815 boundary_functions, 1816 coefficient, 1817 component_mapping); 1818 }, 1819 [&boundary_functions, &dof_to_boundary_mapping, &matrix, &rhs_vector]( 1820 MatrixCreator::internal::AssemblerBoundary :: 1821 CopyData<dim, spacedim, number> const ©_data) { 1822 internal::copy_hp_boundary_mass_matrix_1(copy_data, 1823 boundary_functions, 1824 dof_to_boundary_mapping, 1825 matrix, 1826 rhs_vector); 1827 }, 1828 scratch, 1829 copy_data); 1830 } 1831 1832 1833 1834 template <int dim, int spacedim, typename number> 1835 void create_boundary_mass_matrix(const DoFHandler<dim,spacedim> & dof,const hp::QCollection<dim-1> & q,SparseMatrix<number> & matrix,const std::map<types::boundary_id,const Function<spacedim,number> * > & rhs,Vector<number> & rhs_vector,std::vector<types::global_dof_index> & dof_to_boundary_mapping,const Function<spacedim,number> * const a,std::vector<unsigned int> component_mapping)1836 create_boundary_mass_matrix( 1837 const DoFHandler<dim, spacedim> &dof, 1838 const hp::QCollection<dim - 1> & q, 1839 SparseMatrix<number> & matrix, 1840 const std::map<types::boundary_id, const Function<spacedim, number> *> &rhs, 1841 Vector<number> & rhs_vector, 1842 std::vector<types::global_dof_index> & dof_to_boundary_mapping, 1843 const Function<spacedim, number> *const a, 1844 std::vector<unsigned int> component_mapping) 1845 { 1846 create_boundary_mass_matrix( 1847 hp::StaticMappingQ1<dim, spacedim>::mapping_collection, 1848 dof, 1849 q, 1850 matrix, 1851 rhs, 1852 rhs_vector, 1853 dof_to_boundary_mapping, 1854 a, 1855 component_mapping); 1856 } 1857 1858 1859 1860 template <int dim, int spacedim> 1861 void create_laplace_matrix(const Mapping<dim,spacedim> & mapping,const DoFHandler<dim,spacedim> & dof,const Quadrature<dim> & q,SparseMatrix<double> & matrix,const Function<spacedim> * const coefficient,const AffineConstraints<double> & constraints)1862 create_laplace_matrix(const Mapping<dim, spacedim> & mapping, 1863 const DoFHandler<dim, spacedim> &dof, 1864 const Quadrature<dim> & q, 1865 SparseMatrix<double> & matrix, 1866 const Function<spacedim> *const coefficient, 1867 const AffineConstraints<double> &constraints) 1868 { 1869 Assert(matrix.m() == dof.n_dofs(), 1870 ExcDimensionMismatch(matrix.m(), dof.n_dofs())); 1871 Assert(matrix.n() == dof.n_dofs(), 1872 ExcDimensionMismatch(matrix.n(), dof.n_dofs())); 1873 1874 hp::FECollection<dim, spacedim> fe_collection(dof.get_fe()); 1875 hp::QCollection<dim> q_collection(q); 1876 hp::MappingCollection<dim, spacedim> mapping_collection(mapping); 1877 MatrixCreator::internal::AssemblerData::Scratch<dim, spacedim, double> 1878 assembler_data(fe_collection, 1879 update_gradients | update_JxW_values | 1880 (coefficient != nullptr ? update_quadrature_points : 1881 UpdateFlags(0)), 1882 coefficient, 1883 /*rhs_function=*/nullptr, 1884 q_collection, 1885 mapping_collection); 1886 MatrixCreator::internal::AssemblerData::CopyData<double> copy_data; 1887 copy_data.cell_matrix.reinit( 1888 assembler_data.fe_collection.max_dofs_per_cell(), 1889 assembler_data.fe_collection.max_dofs_per_cell()); 1890 copy_data.cell_rhs.reinit(assembler_data.fe_collection.max_dofs_per_cell()); 1891 copy_data.dof_indices.resize( 1892 assembler_data.fe_collection.max_dofs_per_cell()); 1893 copy_data.constraints = &constraints; 1894 1895 WorkStream::run( 1896 dof.begin_active(), 1897 static_cast<typename DoFHandler<dim, spacedim>::active_cell_iterator>( 1898 dof.end()), 1899 &MatrixCreator::internal::laplace_assembler< 1900 dim, 1901 spacedim, 1902 typename DoFHandler<dim, spacedim>::active_cell_iterator>, 1903 [&matrix](const internal::AssemblerData::CopyData<double> &data) { 1904 MatrixCreator::internal::copy_local_to_global( 1905 data, &matrix, static_cast<Vector<double> *>(nullptr)); 1906 }, 1907 assembler_data, 1908 copy_data); 1909 } 1910 1911 1912 1913 template <int dim, int spacedim> 1914 void create_laplace_matrix(const DoFHandler<dim,spacedim> & dof,const Quadrature<dim> & q,SparseMatrix<double> & matrix,const Function<spacedim> * const coefficient,const AffineConstraints<double> & constraints)1915 create_laplace_matrix(const DoFHandler<dim, spacedim> &dof, 1916 const Quadrature<dim> & q, 1917 SparseMatrix<double> & matrix, 1918 const Function<spacedim> *const coefficient, 1919 const AffineConstraints<double> &constraints) 1920 { 1921 create_laplace_matrix(StaticMappingQ1<dim, spacedim>::mapping, 1922 dof, 1923 q, 1924 matrix, 1925 coefficient, 1926 constraints); 1927 } 1928 1929 1930 1931 template <int dim, int spacedim> 1932 void create_laplace_matrix(const Mapping<dim,spacedim> & mapping,const DoFHandler<dim,spacedim> & dof,const Quadrature<dim> & q,SparseMatrix<double> & matrix,const Function<spacedim> & rhs,Vector<double> & rhs_vector,const Function<spacedim> * const coefficient,const AffineConstraints<double> & constraints)1933 create_laplace_matrix(const Mapping<dim, spacedim> & mapping, 1934 const DoFHandler<dim, spacedim> &dof, 1935 const Quadrature<dim> & q, 1936 SparseMatrix<double> & matrix, 1937 const Function<spacedim> & rhs, 1938 Vector<double> & rhs_vector, 1939 const Function<spacedim> *const coefficient, 1940 const AffineConstraints<double> &constraints) 1941 { 1942 Assert(matrix.m() == dof.n_dofs(), 1943 ExcDimensionMismatch(matrix.m(), dof.n_dofs())); 1944 Assert(matrix.n() == dof.n_dofs(), 1945 ExcDimensionMismatch(matrix.n(), dof.n_dofs())); 1946 1947 hp::FECollection<dim, spacedim> fe_collection(dof.get_fe()); 1948 hp::QCollection<dim> q_collection(q); 1949 hp::MappingCollection<dim, spacedim> mapping_collection(mapping); 1950 MatrixCreator::internal::AssemblerData::Scratch<dim, spacedim, double> 1951 assembler_data(fe_collection, 1952 update_gradients | update_values | update_JxW_values | 1953 update_quadrature_points, 1954 coefficient, 1955 &rhs, 1956 q_collection, 1957 mapping_collection); 1958 MatrixCreator::internal::AssemblerData::CopyData<double> copy_data; 1959 copy_data.cell_matrix.reinit( 1960 assembler_data.fe_collection.max_dofs_per_cell(), 1961 assembler_data.fe_collection.max_dofs_per_cell()); 1962 copy_data.cell_rhs.reinit(assembler_data.fe_collection.max_dofs_per_cell()); 1963 copy_data.dof_indices.resize( 1964 assembler_data.fe_collection.max_dofs_per_cell()); 1965 copy_data.constraints = &constraints; 1966 1967 WorkStream::run( 1968 dof.begin_active(), 1969 static_cast<typename DoFHandler<dim, spacedim>::active_cell_iterator>( 1970 dof.end()), 1971 &MatrixCreator::internal::laplace_assembler< 1972 dim, 1973 spacedim, 1974 typename DoFHandler<dim, spacedim>::active_cell_iterator>, 1975 [&matrix, 1976 &rhs_vector](const internal::AssemblerData::CopyData<double> &data) { 1977 MatrixCreator::internal::copy_local_to_global(data, 1978 &matrix, 1979 &rhs_vector); 1980 }, 1981 assembler_data, 1982 copy_data); 1983 } 1984 1985 1986 1987 template <int dim, int spacedim> 1988 void create_laplace_matrix(const DoFHandler<dim,spacedim> & dof,const Quadrature<dim> & q,SparseMatrix<double> & matrix,const Function<spacedim> & rhs,Vector<double> & rhs_vector,const Function<spacedim> * const coefficient,const AffineConstraints<double> & constraints)1989 create_laplace_matrix(const DoFHandler<dim, spacedim> &dof, 1990 const Quadrature<dim> & q, 1991 SparseMatrix<double> & matrix, 1992 const Function<spacedim> & rhs, 1993 Vector<double> & rhs_vector, 1994 const Function<spacedim> *const coefficient, 1995 const AffineConstraints<double> &constraints) 1996 { 1997 create_laplace_matrix(StaticMappingQ1<dim, spacedim>::mapping, 1998 dof, 1999 q, 2000 matrix, 2001 rhs, 2002 rhs_vector, 2003 coefficient, 2004 constraints); 2005 } 2006 2007 2008 2009 template <int dim, int spacedim> 2010 void create_laplace_matrix(const hp::MappingCollection<dim,spacedim> & mapping,const DoFHandler<dim,spacedim> & dof,const hp::QCollection<dim> & q,SparseMatrix<double> & matrix,const Function<spacedim> * const coefficient,const AffineConstraints<double> & constraints)2011 create_laplace_matrix(const hp::MappingCollection<dim, spacedim> &mapping, 2012 const DoFHandler<dim, spacedim> & dof, 2013 const hp::QCollection<dim> & q, 2014 SparseMatrix<double> & matrix, 2015 const Function<spacedim> *const coefficient, 2016 const AffineConstraints<double> & constraints) 2017 { 2018 Assert(matrix.m() == dof.n_dofs(), 2019 ExcDimensionMismatch(matrix.m(), dof.n_dofs())); 2020 Assert(matrix.n() == dof.n_dofs(), 2021 ExcDimensionMismatch(matrix.n(), dof.n_dofs())); 2022 2023 MatrixCreator::internal::AssemblerData::Scratch<dim, spacedim, double> 2024 assembler_data(dof.get_fe_collection(), 2025 update_gradients | update_JxW_values | 2026 (coefficient != nullptr ? update_quadrature_points : 2027 UpdateFlags(0)), 2028 coefficient, 2029 /*rhs_function=*/nullptr, 2030 q, 2031 mapping); 2032 MatrixCreator::internal::AssemblerData::CopyData<double> copy_data; 2033 copy_data.cell_matrix.reinit( 2034 assembler_data.fe_collection.max_dofs_per_cell(), 2035 assembler_data.fe_collection.max_dofs_per_cell()); 2036 copy_data.cell_rhs.reinit(assembler_data.fe_collection.max_dofs_per_cell()); 2037 copy_data.dof_indices.resize( 2038 assembler_data.fe_collection.max_dofs_per_cell()); 2039 copy_data.constraints = &constraints; 2040 2041 WorkStream::run( 2042 dof.begin_active(), 2043 static_cast<typename DoFHandler<dim, spacedim>::active_cell_iterator>( 2044 dof.end()), 2045 &MatrixCreator::internal::laplace_assembler< 2046 dim, 2047 spacedim, 2048 typename DoFHandler<dim, spacedim>::active_cell_iterator>, 2049 [&matrix](const internal::AssemblerData::CopyData<double> &data) { 2050 MatrixCreator::internal::copy_local_to_global( 2051 data, &matrix, static_cast<Vector<double> *>(nullptr)); 2052 }, 2053 assembler_data, 2054 copy_data); 2055 } 2056 2057 2058 2059 template <int dim, int spacedim> 2060 void create_laplace_matrix(const DoFHandler<dim,spacedim> & dof,const hp::QCollection<dim> & q,SparseMatrix<double> & matrix,const Function<spacedim> * const coefficient,const AffineConstraints<double> & constraints)2061 create_laplace_matrix(const DoFHandler<dim, spacedim> &dof, 2062 const hp::QCollection<dim> & q, 2063 SparseMatrix<double> & matrix, 2064 const Function<spacedim> *const coefficient, 2065 const AffineConstraints<double> &constraints) 2066 { 2067 create_laplace_matrix( 2068 hp::StaticMappingQ1<dim, spacedim>::mapping_collection, 2069 dof, 2070 q, 2071 matrix, 2072 coefficient, 2073 constraints); 2074 } 2075 2076 2077 2078 template <int dim, int spacedim> 2079 void create_laplace_matrix(const hp::MappingCollection<dim,spacedim> & mapping,const DoFHandler<dim,spacedim> & dof,const hp::QCollection<dim> & q,SparseMatrix<double> & matrix,const Function<spacedim> & rhs,Vector<double> & rhs_vector,const Function<spacedim> * const coefficient,const AffineConstraints<double> & constraints)2080 create_laplace_matrix(const hp::MappingCollection<dim, spacedim> &mapping, 2081 const DoFHandler<dim, spacedim> & dof, 2082 const hp::QCollection<dim> & q, 2083 SparseMatrix<double> & matrix, 2084 const Function<spacedim> & rhs, 2085 Vector<double> & rhs_vector, 2086 const Function<spacedim> *const coefficient, 2087 const AffineConstraints<double> & constraints) 2088 { 2089 Assert(matrix.m() == dof.n_dofs(), 2090 ExcDimensionMismatch(matrix.m(), dof.n_dofs())); 2091 Assert(matrix.n() == dof.n_dofs(), 2092 ExcDimensionMismatch(matrix.n(), dof.n_dofs())); 2093 2094 MatrixCreator::internal::AssemblerData::Scratch<dim, spacedim, double> 2095 assembler_data(dof.get_fe_collection(), 2096 update_gradients | update_values | update_JxW_values | 2097 update_quadrature_points, 2098 coefficient, 2099 &rhs, 2100 q, 2101 mapping); 2102 MatrixCreator::internal::AssemblerData::CopyData<double> copy_data; 2103 copy_data.cell_matrix.reinit( 2104 assembler_data.fe_collection.max_dofs_per_cell(), 2105 assembler_data.fe_collection.max_dofs_per_cell()); 2106 copy_data.cell_rhs.reinit(assembler_data.fe_collection.max_dofs_per_cell()); 2107 copy_data.dof_indices.resize( 2108 assembler_data.fe_collection.max_dofs_per_cell()); 2109 copy_data.constraints = &constraints; 2110 2111 WorkStream::run( 2112 dof.begin_active(), 2113 static_cast<typename DoFHandler<dim, spacedim>::active_cell_iterator>( 2114 dof.end()), 2115 &MatrixCreator::internal::laplace_assembler< 2116 dim, 2117 spacedim, 2118 typename DoFHandler<dim, spacedim>::active_cell_iterator>, 2119 [&matrix, 2120 &rhs_vector](const internal::AssemblerData::CopyData<double> &data) { 2121 MatrixCreator::internal::copy_local_to_global(data, 2122 &matrix, 2123 &rhs_vector); 2124 }, 2125 assembler_data, 2126 copy_data); 2127 } 2128 2129 2130 2131 template <int dim, int spacedim> 2132 void create_laplace_matrix(const DoFHandler<dim,spacedim> & dof,const hp::QCollection<dim> & q,SparseMatrix<double> & matrix,const Function<spacedim> & rhs,Vector<double> & rhs_vector,const Function<spacedim> * const coefficient,const AffineConstraints<double> & constraints)2133 create_laplace_matrix(const DoFHandler<dim, spacedim> &dof, 2134 const hp::QCollection<dim> & q, 2135 SparseMatrix<double> & matrix, 2136 const Function<spacedim> & rhs, 2137 Vector<double> & rhs_vector, 2138 const Function<spacedim> *const coefficient, 2139 const AffineConstraints<double> &constraints) 2140 { 2141 create_laplace_matrix( 2142 hp::StaticMappingQ1<dim, spacedim>::mapping_collection, 2143 dof, 2144 q, 2145 matrix, 2146 rhs, 2147 rhs_vector, 2148 coefficient, 2149 constraints); 2150 } 2151 2152 } // namespace MatrixCreator 2153 2154 DEAL_II_NAMESPACE_CLOSE 2155 2156 #endif 2157