1 // --------------------------------------------------------------------- 2 // 3 // Copyright (C) 1999 - 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 #include <deal.II/base/quadrature_lib.h> 17 #include <deal.II/base/table.h> 18 #include <deal.II/base/template_constraints.h> 19 #include <deal.II/base/utilities.h> 20 21 #include <deal.II/distributed/shared_tria.h> 22 #include <deal.II/distributed/tria.h> 23 24 #include <deal.II/dofs/dof_accessor.h> 25 #include <deal.II/dofs/dof_handler.h> 26 #include <deal.II/dofs/dof_tools.h> 27 28 #include <deal.II/fe/fe.h> 29 #include <deal.II/fe/fe_tools.h> 30 #include <deal.II/fe/fe_values.h> 31 32 #include <deal.II/grid/filtered_iterator.h> 33 #include <deal.II/grid/grid_tools.h> 34 #include <deal.II/grid/intergrid_map.h> 35 #include <deal.II/grid/tria.h> 36 #include <deal.II/grid/tria_iterator.h> 37 38 #include <deal.II/hp/dof_handler.h> 39 #include <deal.II/hp/fe_collection.h> 40 #include <deal.II/hp/fe_values.h> 41 #include <deal.II/hp/mapping_collection.h> 42 #include <deal.II/hp/q_collection.h> 43 44 #include <deal.II/lac/affine_constraints.h> 45 #include <deal.II/lac/block_sparsity_pattern.h> 46 #include <deal.II/lac/dynamic_sparsity_pattern.h> 47 #include <deal.II/lac/sparsity_pattern.h> 48 #include <deal.II/lac/trilinos_sparsity_pattern.h> 49 #include <deal.II/lac/vector.h> 50 51 #include <algorithm> 52 #include <numeric> 53 54 DEAL_II_NAMESPACE_OPEN 55 56 57 58 namespace DoFTools 59 { 60 namespace internal 61 { 62 /** 63 * Comparison functor struct to compare two Points and return if one is 64 * "less" than the other one. This can be used to use Point<dim> as a key in 65 * std::map. 66 * 67 * Comparison is done through an artificial downstream direction that 68 * tells directions apart through a factor of 1e-5. Once we got the 69 * direction, we check for its value. In case the distance is exactly zero 70 * (without an epsilon), we might still have the case that two points 71 * combine in a particular way, e.g. the points (1.0, 1.0) and (1.0+1e-5, 72 * 0.0). Thus, compare the points component by component in the second 73 * step. Thus, points need to have identical floating point components to 74 * be considered equal. 75 */ 76 template <int dim, typename Number = double> 77 struct ComparisonHelper 78 { 79 /** 80 * Comparison operator. 81 * 82 * Return true if @p lhs is considered less than @p rhs. 83 */ 84 bool operator ()DoFTools::internal::ComparisonHelper85 operator()(const Point<dim, Number> &lhs, 86 const Point<dim, Number> &rhs) const 87 { 88 double downstream_size = 0; 89 double weight = 1.; 90 for (unsigned int d = 0; d < dim; ++d) 91 { 92 downstream_size += (rhs[d] - lhs[d]) * weight; 93 weight *= 1e-5; 94 } 95 if (downstream_size < 0) 96 return false; 97 else if (downstream_size > 0) 98 return true; 99 else 100 { 101 for (unsigned int d = 0; d < dim; ++d) 102 { 103 if (lhs[d] == rhs[d]) 104 continue; 105 return lhs[d] < rhs[d]; 106 } 107 return false; 108 } 109 } 110 }; 111 112 113 114 // return an array that for each dof on the reference cell 115 // lists the corresponding vector component. 116 // 117 // if an element is non-primitive then we assign to each degree of freedom 118 // the following component: 119 // - if the nonzero components that belong to a shape function are not 120 // selected in the component_mask, then the shape function is assigned 121 // to the first nonzero vector component that corresponds to this 122 // shape function 123 // - otherwise, the shape function is assigned the first component selected 124 // in the component_mask that corresponds to this shape function 125 template <int dim, int spacedim> 126 std::vector<unsigned char> get_local_component_association(const FiniteElement<dim,spacedim> & fe,const ComponentMask & component_mask)127 get_local_component_association(const FiniteElement<dim, spacedim> &fe, 128 const ComponentMask &component_mask) 129 { 130 std::vector<unsigned char> local_component_association( 131 fe.n_dofs_per_cell(), static_cast<unsigned char>(-1)); 132 133 // compute the component each local dof belongs to. 134 // if the shape function is primitive, then this 135 // is simple and we can just associate it with 136 // what system_to_component_index gives us 137 for (unsigned int i = 0; i < fe.n_dofs_per_cell(); ++i) 138 if (fe.is_primitive(i)) 139 local_component_association[i] = 140 fe.system_to_component_index(i).first; 141 else 142 // if the shape function is not primitive, then either use the 143 // component of the first nonzero component corresponding 144 // to this shape function (if the component is not specified 145 // in the component_mask), or the first component of this block 146 // that is listed in the component_mask (if the block this 147 // component corresponds to is indeed specified in the component 148 // mask) 149 { 150 const unsigned int first_comp = 151 fe.get_nonzero_components(i).first_selected_component(); 152 153 if ((fe.get_nonzero_components(i) & component_mask) 154 .n_selected_components(fe.n_components()) == 0) 155 local_component_association[i] = first_comp; 156 else 157 // pick the component selected. we know from the previous 'if' 158 // that within the components that are nonzero for this 159 // shape function there must be at least one for which the 160 // mask is true, so we will for sure run into the break() 161 // at one point 162 for (unsigned int c = first_comp; c < fe.n_components(); ++c) 163 if (component_mask[c] == true) 164 { 165 local_component_association[i] = c; 166 break; 167 } 168 } 169 170 Assert(std::find(local_component_association.begin(), 171 local_component_association.end(), 172 static_cast<unsigned char>(-1)) == 173 local_component_association.end(), 174 ExcInternalError()); 175 176 return local_component_association; 177 } 178 179 180 // this internal function assigns to each dof the respective component 181 // of the vector system. 182 // 183 // the output array dofs_by_component lists for each dof the 184 // corresponding vector component. if the DoFHandler is based on a 185 // parallel distributed triangulation then the output array is index by 186 // dof.locally_owned_dofs().index_within_set(indices[i]) 187 // 188 // if an element is non-primitive then we assign to each degree of 189 // freedom the following component: 190 // - if the nonzero components that belong to a shape function are not 191 // selected in the component_mask, then the shape function is assigned 192 // to the first nonzero vector component that corresponds to this 193 // shape function 194 // - otherwise, the shape function is assigned the first component selected 195 // in the component_mask that corresponds to this shape function 196 template <int dim, int spacedim> 197 void get_component_association(const DoFHandler<dim,spacedim> & dof,const ComponentMask & component_mask,std::vector<unsigned char> & dofs_by_component)198 get_component_association(const DoFHandler<dim, spacedim> &dof, 199 const ComponentMask & component_mask, 200 std::vector<unsigned char> &dofs_by_component) 201 { 202 const dealii::hp::FECollection<dim, spacedim> &fe_collection = 203 dof.get_fe_collection(); 204 Assert(fe_collection.n_components() < 256, ExcNotImplemented()); 205 Assert(dofs_by_component.size() == dof.n_locally_owned_dofs(), 206 ExcDimensionMismatch(dofs_by_component.size(), 207 dof.n_locally_owned_dofs())); 208 209 // next set up a table for the degrees of freedom on each of the 210 // cells (regardless of the fact whether it is listed in the 211 // component_select argument or not) 212 // 213 // for each element 'f' of the FECollection, 214 // local_component_association[f][d] then returns the vector 215 // component that degree of freedom 'd' belongs to 216 std::vector<std::vector<unsigned char>> local_component_association( 217 fe_collection.size()); 218 for (unsigned int f = 0; f < fe_collection.size(); ++f) 219 { 220 const FiniteElement<dim, spacedim> &fe = fe_collection[f]; 221 local_component_association[f] = 222 get_local_component_association(fe, component_mask); 223 } 224 225 // then loop over all cells and do the work 226 std::vector<types::global_dof_index> indices; 227 for (const auto &c : dof.active_cell_iterators()) 228 if (c->is_locally_owned()) 229 { 230 const unsigned int fe_index = c->active_fe_index(); 231 const unsigned int dofs_per_cell = c->get_fe().n_dofs_per_cell(); 232 indices.resize(dofs_per_cell); 233 c->get_dof_indices(indices); 234 for (unsigned int i = 0; i < dofs_per_cell; ++i) 235 if (dof.locally_owned_dofs().is_element(indices[i])) 236 dofs_by_component[dof.locally_owned_dofs().index_within_set( 237 indices[i])] = local_component_association[fe_index][i]; 238 } 239 } 240 241 242 // this is the function corresponding to the one above but working on 243 // blocks instead of components. 244 // 245 // the output array dofs_by_block lists for each dof the corresponding 246 // vector block. if the DoFHandler is based on a parallel distributed 247 // triangulation then the output array is index by 248 // dof.locally_owned_dofs().index_within_set(indices[i]) 249 template <int dim, int spacedim> 250 inline void get_block_association(const DoFHandler<dim,spacedim> & dof,std::vector<unsigned char> & dofs_by_block)251 get_block_association(const DoFHandler<dim, spacedim> &dof, 252 std::vector<unsigned char> & dofs_by_block) 253 { 254 const dealii::hp::FECollection<dim, spacedim> &fe_collection = 255 dof.get_fe_collection(); 256 Assert(fe_collection.n_components() < 256, ExcNotImplemented()); 257 Assert(dofs_by_block.size() == dof.n_locally_owned_dofs(), 258 ExcDimensionMismatch(dofs_by_block.size(), 259 dof.n_locally_owned_dofs())); 260 261 // next set up a table for the degrees of freedom on each of the 262 // cells (regardless of the fact whether it is listed in the 263 // component_select argument or not) 264 // 265 // for each element 'f' of the FECollection, 266 // local_block_association[f][d] then returns the vector block that 267 // degree of freedom 'd' belongs to 268 std::vector<std::vector<unsigned char>> local_block_association( 269 fe_collection.size()); 270 for (unsigned int f = 0; f < fe_collection.size(); ++f) 271 { 272 const FiniteElement<dim, spacedim> &fe = fe_collection[f]; 273 local_block_association[f].resize(fe.n_dofs_per_cell(), 274 static_cast<unsigned char>(-1)); 275 for (unsigned int i = 0; i < fe.n_dofs_per_cell(); ++i) 276 local_block_association[f][i] = fe.system_to_block_index(i).first; 277 278 Assert(std::find(local_block_association[f].begin(), 279 local_block_association[f].end(), 280 static_cast<unsigned char>(-1)) == 281 local_block_association[f].end(), 282 ExcInternalError()); 283 } 284 285 // then loop over all cells and do the work 286 std::vector<types::global_dof_index> indices; 287 for (const auto &cell : dof.active_cell_iterators()) 288 if (cell->is_locally_owned()) 289 { 290 const unsigned int fe_index = cell->active_fe_index(); 291 const unsigned int dofs_per_cell = cell->get_fe().n_dofs_per_cell(); 292 indices.resize(dofs_per_cell); 293 cell->get_dof_indices(indices); 294 for (unsigned int i = 0; i < dofs_per_cell; ++i) 295 if (dof.locally_owned_dofs().is_element(indices[i])) 296 dofs_by_block[dof.locally_owned_dofs().index_within_set( 297 indices[i])] = local_block_association[fe_index][i]; 298 } 299 } 300 } // namespace internal 301 302 303 304 template <int dim, int spacedim, typename Number> 305 void distribute_cell_to_dof_vector(const DoFHandler<dim,spacedim> & dof_handler,const Vector<Number> & cell_data,Vector<double> & dof_data,const unsigned int component)306 distribute_cell_to_dof_vector(const DoFHandler<dim, spacedim> &dof_handler, 307 const Vector<Number> & cell_data, 308 Vector<double> & dof_data, 309 const unsigned int component) 310 { 311 const Triangulation<dim, spacedim> &tria = dof_handler.get_triangulation(); 312 (void)tria; 313 314 AssertDimension(cell_data.size(), tria.n_active_cells()); 315 AssertDimension(dof_data.size(), dof_handler.n_dofs()); 316 const auto &fe_collection = dof_handler.get_fe_collection(); 317 AssertIndexRange(component, fe_collection.n_components()); 318 for (unsigned int i = 0; i < fe_collection.size(); ++i) 319 { 320 Assert(fe_collection[i].is_primitive() == true, 321 typename FiniteElement<dim>::ExcFENotPrimitive()); 322 } 323 324 // store a flag whether we should care about different components. this 325 // is just a simplification, we could ask for this at every single 326 // place equally well 327 const bool consider_components = 328 (dof_handler.get_fe_collection().n_components() != 1); 329 330 // zero out the components that we will touch 331 if (consider_components == false) 332 dof_data = 0; 333 else 334 { 335 std::vector<unsigned char> component_dofs( 336 dof_handler.n_locally_owned_dofs()); 337 internal::get_component_association( 338 dof_handler, 339 fe_collection.component_mask(FEValuesExtractors::Scalar(component)), 340 component_dofs); 341 342 for (unsigned int i = 0; i < dof_data.size(); ++i) 343 if (component_dofs[i] == static_cast<unsigned char>(component)) 344 dof_data(i) = 0; 345 } 346 347 // count how often we have added a value in the sum for each dof 348 std::vector<unsigned char> touch_count(dof_handler.n_dofs(), 0); 349 350 typename DoFHandler<dim, spacedim>::active_cell_iterator 351 cell = dof_handler.begin_active(), 352 endc = dof_handler.end(); 353 std::vector<types::global_dof_index> dof_indices; 354 dof_indices.reserve(fe_collection.max_dofs_per_cell()); 355 356 for (unsigned int present_cell = 0; cell != endc; ++cell, ++present_cell) 357 { 358 const unsigned int dofs_per_cell = cell->get_fe().n_dofs_per_cell(); 359 dof_indices.resize(dofs_per_cell); 360 cell->get_dof_indices(dof_indices); 361 362 for (unsigned int i = 0; i < dofs_per_cell; ++i) 363 // consider this dof only if it is the right component. if there 364 // is only one component, short cut the test 365 if (!consider_components || 366 (cell->get_fe().system_to_component_index(i).first == component)) 367 { 368 // sum up contribution of the present_cell to this dof 369 dof_data(dof_indices[i]) += cell_data(present_cell); 370 371 // note that we added another summand 372 ++touch_count[dof_indices[i]]; 373 } 374 } 375 376 // compute the mean value on all the dofs by dividing with the number 377 // of summands. 378 for (types::global_dof_index i = 0; i < dof_handler.n_dofs(); ++i) 379 { 380 // assert that each dof was used at least once. this needs not be 381 // the case if the vector has more than one component 382 Assert(consider_components || (touch_count[i] != 0), 383 ExcInternalError()); 384 if (touch_count[i] != 0) 385 dof_data(i) /= touch_count[i]; 386 } 387 } 388 389 390 391 template <int dim, int spacedim> 392 void extract_dofs(const DoFHandler<dim,spacedim> & dof,const ComponentMask & component_mask,std::vector<bool> & selected_dofs)393 extract_dofs(const DoFHandler<dim, spacedim> &dof, 394 const ComponentMask & component_mask, 395 std::vector<bool> & selected_dofs) 396 { 397 Assert(component_mask.represents_n_components( 398 dof.get_fe_collection().n_components()), 399 ExcMessage( 400 "The given component mask is not sized correctly to represent the " 401 "components of the given finite element.")); 402 Assert(selected_dofs.size() == dof.n_locally_owned_dofs(), 403 ExcDimensionMismatch(selected_dofs.size(), 404 dof.n_locally_owned_dofs())); 405 406 // two special cases: no component is selected, and all components are 407 // selected; both rather stupid, but easy to catch 408 if (component_mask.n_selected_components( 409 dof.get_fe_collection().n_components()) == 0) 410 { 411 std::fill_n(selected_dofs.begin(), dof.n_locally_owned_dofs(), false); 412 return; 413 } 414 else if (component_mask.n_selected_components( 415 dof.get_fe_collection().n_components()) == 416 dof.get_fe_collection().n_components()) 417 { 418 std::fill_n(selected_dofs.begin(), dof.n_locally_owned_dofs(), true); 419 return; 420 } 421 422 423 // preset all values by false 424 std::fill_n(selected_dofs.begin(), dof.n_locally_owned_dofs(), false); 425 426 // get the component association of each DoF and then select the ones 427 // that match the given set of components 428 std::vector<unsigned char> dofs_by_component(dof.n_locally_owned_dofs()); 429 internal::get_component_association(dof, component_mask, dofs_by_component); 430 431 for (types::global_dof_index i = 0; i < dof.n_locally_owned_dofs(); ++i) 432 if (component_mask[dofs_by_component[i]] == true) 433 selected_dofs[i] = true; 434 } 435 436 437 438 template <int dim, int spacedim> 439 IndexSet extract_dofs(const DoFHandler<dim,spacedim> & dof,const ComponentMask & component_mask)440 extract_dofs(const DoFHandler<dim, spacedim> &dof, 441 const ComponentMask & component_mask) 442 { 443 Assert(component_mask.represents_n_components( 444 dof.get_fe_collection().n_components()), 445 ExcMessage( 446 "The given component mask is not sized correctly to represent the " 447 "components of the given finite element.")); 448 449 // Two special cases: no component is selected, and all components are 450 // selected; both rather stupid, but easy to catch 451 if (component_mask.n_selected_components( 452 dof.get_fe_collection().n_components()) == 0) 453 return IndexSet(dof.n_dofs()); 454 else if (component_mask.n_selected_components( 455 dof.get_fe_collection().n_components()) == 456 dof.get_fe_collection().n_components()) 457 return dof.locally_owned_dofs(); 458 459 // get the component association of each DoF and then select the ones 460 // that match the given set of components 461 std::vector<unsigned char> dofs_by_component(dof.n_locally_owned_dofs()); 462 internal::get_component_association(dof, component_mask, dofs_by_component); 463 464 // fill the selected components in a vector 465 std::vector<types::global_dof_index> selected_dofs; 466 selected_dofs.reserve(dof.n_locally_owned_dofs()); 467 for (types::global_dof_index i = 0; i < dofs_by_component.size(); ++i) 468 if (component_mask[dofs_by_component[i]] == true) 469 selected_dofs.push_back(dof.locally_owned_dofs().nth_index_in_set(i)); 470 471 // fill vector of indices to return argument 472 IndexSet result(dof.n_dofs()); 473 result.add_indices(selected_dofs.begin(), selected_dofs.end()); 474 return result; 475 } 476 477 478 479 template <int dim, int spacedim> 480 void extract_dofs(const DoFHandler<dim,spacedim> & dof,const BlockMask & block_mask,std::vector<bool> & selected_dofs)481 extract_dofs(const DoFHandler<dim, spacedim> &dof, 482 const BlockMask & block_mask, 483 std::vector<bool> & selected_dofs) 484 { 485 // simply forward to the function that works based on a component mask 486 extract_dofs<dim, spacedim>( 487 dof, dof.get_fe_collection().component_mask(block_mask), selected_dofs); 488 } 489 490 491 492 template <int dim, int spacedim> 493 IndexSet extract_dofs(const DoFHandler<dim,spacedim> & dof,const BlockMask & block_mask)494 extract_dofs(const DoFHandler<dim, spacedim> &dof, 495 const BlockMask & block_mask) 496 { 497 // simply forward to the function that works based on a component mask 498 return extract_dofs<dim, spacedim>( 499 dof, dof.get_fe_collection().component_mask(block_mask)); 500 } 501 502 503 504 template <int dim, int spacedim> 505 std::vector<IndexSet> locally_owned_dofs_per_component(const DoFHandler<dim,spacedim> & dof,const ComponentMask & component_mask)506 locally_owned_dofs_per_component(const DoFHandler<dim, spacedim> &dof, 507 const ComponentMask &component_mask) 508 { 509 const auto n_comps = dof.get_fe_collection().n_components(); 510 Assert(component_mask.represents_n_components(n_comps), 511 ExcMessage( 512 "The given component mask is not sized correctly to represent the " 513 "components of the given finite element.")); 514 515 const auto &locally_owned_dofs = dof.locally_owned_dofs(); 516 517 // get the component association of each DoF and then select the ones 518 // that match the given set of components 519 std::vector<unsigned char> dofs_by_component(dof.n_locally_owned_dofs()); 520 internal::get_component_association(dof, component_mask, dofs_by_component); 521 522 std::vector<IndexSet> index_per_comp(n_comps, IndexSet(dof.n_dofs())); 523 524 for (types::global_dof_index i = 0; i < dof.n_locally_owned_dofs(); ++i) 525 { 526 const auto &comp_i = dofs_by_component[i]; 527 if (component_mask[comp_i]) 528 index_per_comp[comp_i].add_index( 529 locally_owned_dofs.nth_index_in_set(i)); 530 } 531 for (auto &c : index_per_comp) 532 c.compress(); 533 return index_per_comp; 534 } 535 536 537 538 template <int dim, int spacedim> 539 void extract_level_dofs(const unsigned int level,const DoFHandler<dim,spacedim> & dof,const ComponentMask & component_mask,std::vector<bool> & selected_dofs)540 extract_level_dofs(const unsigned int level, 541 const DoFHandler<dim, spacedim> &dof, 542 const ComponentMask & component_mask, 543 std::vector<bool> & selected_dofs) 544 { 545 const FiniteElement<dim, spacedim> &fe = dof.get_fe(); 546 547 Assert(component_mask.represents_n_components( 548 dof.get_fe_collection().n_components()), 549 ExcMessage( 550 "The given component mask is not sized correctly to represent the " 551 "components of the given finite element.")); 552 Assert(selected_dofs.size() == dof.n_dofs(level), 553 ExcDimensionMismatch(selected_dofs.size(), dof.n_dofs(level))); 554 555 // two special cases: no component is selected, and all components are 556 // selected, both rather stupid, but easy to catch 557 if (component_mask.n_selected_components( 558 dof.get_fe_collection().n_components()) == 0) 559 { 560 std::fill_n(selected_dofs.begin(), dof.n_dofs(level), false); 561 return; 562 } 563 else if (component_mask.n_selected_components( 564 dof.get_fe_collection().n_components()) == 565 dof.get_fe_collection().n_components()) 566 { 567 std::fill_n(selected_dofs.begin(), dof.n_dofs(level), true); 568 return; 569 } 570 571 // preset all values by false 572 std::fill_n(selected_dofs.begin(), dof.n_dofs(level), false); 573 574 // next set up a table for the degrees of freedom on each of the cells 575 // whether it is something interesting or not 576 std::vector<unsigned char> local_component_asssociation = 577 internal::get_local_component_association(fe, component_mask); 578 std::vector<bool> local_selected_dofs(fe.n_dofs_per_cell()); 579 for (unsigned int i = 0; i < fe.n_dofs_per_cell(); ++i) 580 local_selected_dofs[i] = component_mask[local_component_asssociation[i]]; 581 582 // then loop over all cells and do work 583 std::vector<types::global_dof_index> indices(fe.n_dofs_per_cell()); 584 typename DoFHandler<dim, spacedim>::level_cell_iterator c; 585 for (c = dof.begin(level); c != dof.end(level); ++c) 586 { 587 c->get_mg_dof_indices(indices); 588 for (unsigned int i = 0; i < fe.n_dofs_per_cell(); ++i) 589 selected_dofs[indices[i]] = local_selected_dofs[i]; 590 } 591 } 592 593 594 595 template <int dim, int spacedim> 596 void extract_level_dofs(const unsigned int level,const DoFHandler<dim,spacedim> & dof,const BlockMask & block_mask,std::vector<bool> & selected_dofs)597 extract_level_dofs(const unsigned int level, 598 const DoFHandler<dim, spacedim> &dof, 599 const BlockMask & block_mask, 600 std::vector<bool> & selected_dofs) 601 { 602 // simply defer to the other extract_level_dofs() function 603 extract_level_dofs(level, 604 dof, 605 dof.get_fe().component_mask(block_mask), 606 selected_dofs); 607 } 608 609 610 611 template <int dim, int spacedim> 612 void extract_boundary_dofs(const DoFHandler<dim,spacedim> & dof_handler,const ComponentMask & component_mask,std::vector<bool> & selected_dofs,const std::set<types::boundary_id> & boundary_ids)613 extract_boundary_dofs(const DoFHandler<dim, spacedim> & dof_handler, 614 const ComponentMask & component_mask, 615 std::vector<bool> & selected_dofs, 616 const std::set<types::boundary_id> &boundary_ids) 617 { 618 Assert((dynamic_cast< 619 const parallel::distributed::Triangulation<dim, spacedim> *>( 620 &dof_handler.get_triangulation()) == nullptr), 621 ExcMessage( 622 "This function can not be used with distributed triangulations. " 623 "See the documentation for more information.")); 624 625 IndexSet indices; 626 extract_boundary_dofs(dof_handler, component_mask, indices, boundary_ids); 627 628 // clear and reset array by default values 629 selected_dofs.clear(); 630 selected_dofs.resize(dof_handler.n_dofs(), false); 631 632 // then convert the values computed above to the binary vector 633 indices.fill_binary_vector(selected_dofs); 634 } 635 636 637 template <int dim, int spacedim> 638 void extract_boundary_dofs(const DoFHandler<dim,spacedim> & dof_handler,const ComponentMask & component_mask,IndexSet & selected_dofs,const std::set<types::boundary_id> & boundary_ids)639 extract_boundary_dofs(const DoFHandler<dim, spacedim> & dof_handler, 640 const ComponentMask & component_mask, 641 IndexSet & selected_dofs, 642 const std::set<types::boundary_id> &boundary_ids) 643 { 644 Assert(component_mask.represents_n_components( 645 dof_handler.get_fe_collection().n_components()), 646 ExcMessage("Component mask has invalid size.")); 647 Assert(boundary_ids.find(numbers::internal_face_boundary_id) == 648 boundary_ids.end(), 649 ExcInvalidBoundaryIndicator()); 650 651 // first reset output argument 652 selected_dofs.clear(); 653 selected_dofs.set_size(dof_handler.n_dofs()); 654 655 // let's see whether we have to check for certain boundary indicators 656 // or whether we can accept all 657 const bool check_boundary_id = (boundary_ids.size() != 0); 658 659 // also see whether we have to check whether a certain vector component 660 // is selected, or all 661 const bool check_vector_component = 662 ((component_mask.represents_the_all_selected_mask() == false) || 663 (component_mask.n_selected_components( 664 dof_handler.get_fe_collection().n_components()) != 665 dof_handler.get_fe_collection().n_components())); 666 667 std::vector<types::global_dof_index> face_dof_indices; 668 face_dof_indices.reserve( 669 dof_handler.get_fe_collection().max_dofs_per_face()); 670 671 // now loop over all cells and check whether their faces are at the 672 // boundary. note that we need not take special care of single lines 673 // being at the boundary (using @p{cell->has_boundary_lines}), since we 674 // do not support boundaries of dimension dim-2, and so every isolated 675 // boundary line is also part of a boundary face which we will be 676 // visiting sooner or later 677 for (const auto &cell : dof_handler.active_cell_iterators()) 678 679 // only work on cells that are either locally owned or at least ghost 680 // cells 681 if (cell->is_artificial() == false) 682 for (const unsigned int face : cell->face_indices()) 683 if (cell->at_boundary(face)) 684 if (!check_boundary_id || 685 (boundary_ids.find(cell->face(face)->boundary_id()) != 686 boundary_ids.end())) 687 { 688 const FiniteElement<dim, spacedim> &fe = cell->get_fe(); 689 690 const unsigned int dofs_per_face = fe.n_dofs_per_face(); 691 face_dof_indices.resize(dofs_per_face); 692 cell->face(face)->get_dof_indices(face_dof_indices, 693 cell->active_fe_index()); 694 695 for (unsigned int i = 0; i < fe.n_dofs_per_face(); ++i) 696 if (!check_vector_component) 697 selected_dofs.add_index(face_dof_indices[i]); 698 else 699 // check for component is required. somewhat tricky as 700 // usual for the case that the shape function is 701 // non-primitive, but use usual convention (see docs) 702 { 703 // first get at the cell-global number of a face dof, 704 // to ask the fe certain questions 705 const unsigned int cell_index = 706 (dim == 1 ? 707 i : 708 (dim == 2 ? 709 (i < 2 * fe.n_dofs_per_vertex() ? 710 i : 711 i + 2 * fe.n_dofs_per_vertex()) : 712 (dim == 3 ? (i < 4 * fe.n_dofs_per_vertex() ? 713 i : 714 (i < 4 * fe.n_dofs_per_vertex() + 715 4 * fe.n_dofs_per_line() ? 716 i + 4 * fe.n_dofs_per_vertex() : 717 i + 4 * fe.n_dofs_per_vertex() + 718 8 * fe.n_dofs_per_line())) : 719 numbers::invalid_unsigned_int))); 720 if (fe.is_primitive(cell_index)) 721 { 722 if (component_mask 723 [fe.face_system_to_component_index(i).first] == 724 true) 725 selected_dofs.add_index(face_dof_indices[i]); 726 } 727 else // not primitive 728 { 729 const unsigned int first_nonzero_comp = 730 fe.get_nonzero_components(cell_index) 731 .first_selected_component(); 732 Assert(first_nonzero_comp < fe.n_components(), 733 ExcInternalError()); 734 735 if (component_mask[first_nonzero_comp] == true) 736 selected_dofs.add_index(face_dof_indices[i]); 737 } 738 } 739 } 740 } 741 742 743 744 template <int dim, int spacedim> 745 void extract_dofs_with_support_on_boundary(const DoFHandler<dim,spacedim> & dof_handler,const ComponentMask & component_mask,std::vector<bool> & selected_dofs,const std::set<types::boundary_id> & boundary_ids)746 extract_dofs_with_support_on_boundary( 747 const DoFHandler<dim, spacedim> & dof_handler, 748 const ComponentMask & component_mask, 749 std::vector<bool> & selected_dofs, 750 const std::set<types::boundary_id> &boundary_ids) 751 { 752 Assert(component_mask.represents_n_components( 753 dof_handler.get_fe_collection().n_components()), 754 ExcMessage("This component mask has the wrong size.")); 755 Assert(boundary_ids.find(numbers::internal_face_boundary_id) == 756 boundary_ids.end(), 757 ExcInvalidBoundaryIndicator()); 758 759 // let's see whether we have to check for certain boundary indicators 760 // or whether we can accept all 761 const bool check_boundary_id = (boundary_ids.size() != 0); 762 763 // also see whether we have to check whether a certain vector component 764 // is selected, or all 765 const bool check_vector_component = 766 (component_mask.represents_the_all_selected_mask() == false); 767 768 // clear and reset array by default values 769 selected_dofs.clear(); 770 selected_dofs.resize(dof_handler.n_dofs(), false); 771 std::vector<types::global_dof_index> cell_dof_indices; 772 cell_dof_indices.reserve( 773 dof_handler.get_fe_collection().max_dofs_per_cell()); 774 775 // now loop over all cells and check whether their faces are at the 776 // boundary. note that we need not take special care of single lines 777 // being at the boundary (using @p{cell->has_boundary_lines}), since we 778 // do not support boundaries of dimension dim-2, and so every isolated 779 // boundary line is also part of a boundary face which we will be 780 // visiting sooner or later 781 for (const auto &cell : dof_handler.active_cell_iterators()) 782 for (const unsigned int face : cell->face_indices()) 783 if (cell->at_boundary(face)) 784 if (!check_boundary_id || 785 (boundary_ids.find(cell->face(face)->boundary_id()) != 786 boundary_ids.end())) 787 { 788 const FiniteElement<dim, spacedim> &fe = cell->get_fe(); 789 790 const unsigned int dofs_per_cell = fe.n_dofs_per_cell(); 791 cell_dof_indices.resize(dofs_per_cell); 792 cell->get_dof_indices(cell_dof_indices); 793 794 for (unsigned int i = 0; i < fe.n_dofs_per_cell(); ++i) 795 if (fe.has_support_on_face(i, face)) 796 { 797 if (!check_vector_component) 798 selected_dofs[cell_dof_indices[i]] = true; 799 else 800 // check for component is required. somewhat tricky 801 // as usual for the case that the shape function is 802 // non-primitive, but use usual convention (see docs) 803 { 804 if (fe.is_primitive(i)) 805 selected_dofs[cell_dof_indices[i]] = 806 (component_mask[fe.system_to_component_index(i) 807 .first] == true); 808 else // not primitive 809 { 810 const unsigned int first_nonzero_comp = 811 fe.get_nonzero_components(i) 812 .first_selected_component(); 813 Assert(first_nonzero_comp < fe.n_components(), 814 ExcInternalError()); 815 816 selected_dofs[cell_dof_indices[i]] = 817 (component_mask[first_nonzero_comp] == true); 818 } 819 } 820 } 821 } 822 } 823 824 825 826 template <int dim, int spacedim, typename number> 827 IndexSet extract_dofs_with_support_contained_within(const DoFHandler<dim,spacedim> & dof_handler,const std::function<bool (const typename DoFHandler<dim,spacedim>::active_cell_iterator &)> & predicate,const AffineConstraints<number> & cm)828 extract_dofs_with_support_contained_within( 829 const DoFHandler<dim, spacedim> &dof_handler, 830 const std::function< 831 bool(const typename DoFHandler<dim, spacedim>::active_cell_iterator &)> 832 & predicate, 833 const AffineConstraints<number> &cm) 834 { 835 const std::function<bool( 836 const typename DoFHandler<dim, spacedim>::active_cell_iterator &)> 837 predicate_local = 838 [=]( 839 const typename DoFHandler<dim, spacedim>::active_cell_iterator &cell) 840 -> bool { return cell->is_locally_owned() && predicate(cell); }; 841 842 std::vector<types::global_dof_index> local_dof_indices; 843 local_dof_indices.reserve( 844 dof_handler.get_fe_collection().max_dofs_per_cell()); 845 846 // Get all the dofs that live on the subdomain: 847 std::set<types::global_dof_index> predicate_dofs; 848 849 typename DoFHandler<dim, spacedim>::active_cell_iterator 850 cell = dof_handler.begin_active(), 851 endc = dof_handler.end(); 852 for (; cell != endc; ++cell) 853 if (!cell->is_artificial() && predicate(cell)) 854 { 855 local_dof_indices.resize(cell->get_fe().n_dofs_per_cell()); 856 cell->get_dof_indices(local_dof_indices); 857 predicate_dofs.insert(local_dof_indices.begin(), 858 local_dof_indices.end()); 859 } 860 861 // Get halo layer and accumulate its DoFs 862 std::set<types::global_dof_index> dofs_with_support_on_halo_cells; 863 864 const std::vector<typename DoFHandler<dim, spacedim>::active_cell_iterator> 865 halo_cells = 866 GridTools::compute_active_cell_halo_layer(dof_handler, predicate_local); 867 for (typename std::vector<typename DoFHandler<dim, spacedim>:: 868 active_cell_iterator>::const_iterator it = 869 halo_cells.begin(); 870 it != halo_cells.end(); 871 ++it) 872 { 873 // skip halo cells that satisfy the given predicate and thereby 874 // shall be a part of the index set on another MPI core. 875 // Those could only be ghost cells with p::d::Tria. 876 if (predicate(*it)) 877 { 878 Assert((*it)->is_ghost(), ExcInternalError()); 879 continue; 880 } 881 882 const unsigned int dofs_per_cell = (*it)->get_fe().n_dofs_per_cell(); 883 local_dof_indices.resize(dofs_per_cell); 884 (*it)->get_dof_indices(local_dof_indices); 885 dofs_with_support_on_halo_cells.insert(local_dof_indices.begin(), 886 local_dof_indices.end()); 887 } 888 889 // A complication coming from the fact that DoFs living on locally 890 // owned cells which satisfy predicate may participate in constraints for 891 // DoFs outside of this region. 892 if (cm.n_constraints() > 0) 893 { 894 // Remove DoFs that are part of constraints for DoFs outside 895 // of predicate. Since we will subtract halo_dofs from predicate_dofs, 896 // achieve this by extending halo_dofs with DoFs to which 897 // halo_dofs are constrained. 898 std::set<types::global_dof_index> extra_halo; 899 for (const auto dof : dofs_with_support_on_halo_cells) 900 // if halo DoF is constrained, add all DoFs to which it's constrained 901 // because after resolving constraints, the support of the DoFs that 902 // constrain the current DoF will extend to the halo cells. 903 if (const auto *line_ptr = cm.get_constraint_entries(dof)) 904 { 905 const unsigned int line_size = line_ptr->size(); 906 for (unsigned int j = 0; j < line_size; ++j) 907 extra_halo.insert((*line_ptr)[j].first); 908 } 909 910 dofs_with_support_on_halo_cells.insert(extra_halo.begin(), 911 extra_halo.end()); 912 } 913 914 // Rework our std::set's into IndexSet and subtract halo DoFs from the 915 // predicate_dofs: 916 IndexSet support_set(dof_handler.n_dofs()); 917 support_set.add_indices(predicate_dofs.begin(), predicate_dofs.end()); 918 support_set.compress(); 919 920 IndexSet halo_set(dof_handler.n_dofs()); 921 halo_set.add_indices(dofs_with_support_on_halo_cells.begin(), 922 dofs_with_support_on_halo_cells.end()); 923 halo_set.compress(); 924 925 support_set.subtract_set(halo_set); 926 927 // we intentionally do not want to limit the output to locally owned DoFs. 928 return support_set; 929 } 930 931 932 933 namespace internal 934 { 935 namespace 936 { 937 template <int spacedim> 938 IndexSet extract_hanging_node_dofs(const dealii::DoFHandler<1,spacedim> & dof_handler)939 extract_hanging_node_dofs( 940 const dealii::DoFHandler<1, spacedim> &dof_handler) 941 { 942 // there are no hanging nodes in 1d 943 return IndexSet(dof_handler.n_dofs()); 944 } 945 946 947 template <int spacedim> 948 IndexSet extract_hanging_node_dofs(const dealii::DoFHandler<2,spacedim> & dof_handler)949 extract_hanging_node_dofs( 950 const dealii::DoFHandler<2, spacedim> &dof_handler) 951 { 952 const unsigned int dim = 2; 953 954 IndexSet selected_dofs(dof_handler.n_dofs()); 955 956 const FiniteElement<dim, spacedim> &fe = dof_handler.get_fe(); 957 958 // this function is similar to the make_sparsity_pattern function, 959 // see there for more information 960 for (const auto &cell : dof_handler.active_cell_iterators()) 961 if (!cell->is_artificial()) 962 { 963 for (const unsigned int face : cell->face_indices()) 964 if (cell->face(face)->has_children()) 965 { 966 const typename dealii::DoFHandler<dim, 967 spacedim>::line_iterator 968 line = cell->face(face); 969 970 for (unsigned int dof = 0; dof != fe.n_dofs_per_vertex(); 971 ++dof) 972 selected_dofs.add_index( 973 line->child(0)->vertex_dof_index(1, dof)); 974 975 for (unsigned int child = 0; child < 2; ++child) 976 { 977 if (cell->neighbor_child_on_subface(face, child) 978 ->is_artificial()) 979 continue; 980 for (unsigned int dof = 0; dof != fe.n_dofs_per_line(); 981 ++dof) 982 selected_dofs.add_index( 983 line->child(child)->dof_index(dof)); 984 } 985 } 986 } 987 988 selected_dofs.compress(); 989 return selected_dofs; 990 } 991 992 993 template <int spacedim> 994 IndexSet extract_hanging_node_dofs(const dealii::DoFHandler<3,spacedim> & dof_handler)995 extract_hanging_node_dofs( 996 const dealii::DoFHandler<3, spacedim> &dof_handler) 997 { 998 const unsigned int dim = 3; 999 1000 IndexSet selected_dofs(dof_handler.n_dofs()); 1001 IndexSet unconstrained_dofs(dof_handler.n_dofs()); 1002 1003 const FiniteElement<dim, spacedim> &fe = dof_handler.get_fe(); 1004 1005 for (const auto &cell : dof_handler.active_cell_iterators()) 1006 if (!cell->is_artificial()) 1007 for (auto f : cell->face_indices()) 1008 { 1009 const typename dealii::DoFHandler<dim, spacedim>::face_iterator 1010 face = cell->face(f); 1011 if (cell->face(f)->has_children()) 1012 { 1013 for (unsigned int child = 0; child < 4; ++child) 1014 if (!cell->neighbor_child_on_subface(f, child) 1015 ->is_artificial()) 1016 { 1017 // simply take all DoFs that live on this subface 1018 std::vector<types::global_dof_index> ldi( 1019 fe.n_dofs_per_face()); 1020 face->child(child)->get_dof_indices(ldi); 1021 selected_dofs.add_indices(ldi.begin(), ldi.end()); 1022 } 1023 1024 // and subtract (in the end) all the indices which a shared 1025 // between this face and its subfaces 1026 for (unsigned int vertex = 0; vertex < 4; ++vertex) 1027 for (unsigned int dof = 0; dof != fe.n_dofs_per_vertex(); 1028 ++dof) 1029 unconstrained_dofs.add_index( 1030 face->vertex_dof_index(vertex, dof)); 1031 } 1032 } 1033 selected_dofs.subtract_set(unconstrained_dofs); 1034 return selected_dofs; 1035 } 1036 } // namespace 1037 } // namespace internal 1038 1039 1040 1041 template <int dim, int spacedim> 1042 void extract_hanging_node_dofs(const DoFHandler<dim,spacedim> & dof_handler,std::vector<bool> & selected_dofs)1043 extract_hanging_node_dofs(const DoFHandler<dim, spacedim> &dof_handler, 1044 std::vector<bool> & selected_dofs) 1045 { 1046 const IndexSet selected_dofs_as_index_set = 1047 extract_hanging_node_dofs(dof_handler); 1048 Assert(selected_dofs.size() == dof_handler.n_dofs(), 1049 ExcDimensionMismatch(selected_dofs.size(), dof_handler.n_dofs())); 1050 // preset all values by false 1051 std::fill(selected_dofs.begin(), selected_dofs.end(), false); 1052 for (const auto index : selected_dofs_as_index_set) 1053 selected_dofs[index] = true; 1054 } 1055 1056 1057 1058 template <int dim, int spacedim> 1059 IndexSet extract_hanging_node_dofs(const DoFHandler<dim,spacedim> & dof_handler)1060 extract_hanging_node_dofs(const DoFHandler<dim, spacedim> &dof_handler) 1061 { 1062 return internal::extract_hanging_node_dofs(dof_handler); 1063 } 1064 1065 1066 1067 template <int dim, int spacedim> 1068 void extract_subdomain_dofs(const DoFHandler<dim,spacedim> & dof_handler,const types::subdomain_id subdomain_id,std::vector<bool> & selected_dofs)1069 extract_subdomain_dofs(const DoFHandler<dim, spacedim> &dof_handler, 1070 const types::subdomain_id subdomain_id, 1071 std::vector<bool> & selected_dofs) 1072 { 1073 Assert(selected_dofs.size() == dof_handler.n_dofs(), 1074 ExcDimensionMismatch(selected_dofs.size(), dof_handler.n_dofs())); 1075 1076 // preset all values by false 1077 std::fill_n(selected_dofs.begin(), dof_handler.n_dofs(), false); 1078 1079 std::vector<types::global_dof_index> local_dof_indices; 1080 local_dof_indices.reserve( 1081 dof_handler.get_fe_collection().max_dofs_per_cell()); 1082 1083 // this function is similar to the make_sparsity_pattern function, see 1084 // there for more information 1085 typename DoFHandler<dim, spacedim>::active_cell_iterator 1086 cell = dof_handler.begin_active(), 1087 endc = dof_handler.end(); 1088 for (; cell != endc; ++cell) 1089 if (cell->subdomain_id() == subdomain_id) 1090 { 1091 const unsigned int dofs_per_cell = cell->get_fe().n_dofs_per_cell(); 1092 local_dof_indices.resize(dofs_per_cell); 1093 cell->get_dof_indices(local_dof_indices); 1094 for (unsigned int i = 0; i < dofs_per_cell; ++i) 1095 selected_dofs[local_dof_indices[i]] = true; 1096 } 1097 } 1098 1099 1100 1101 template <int dim, int spacedim> 1102 void extract_locally_active_dofs(const DoFHandler<dim,spacedim> & dof_handler,IndexSet & dof_set)1103 extract_locally_active_dofs(const DoFHandler<dim, spacedim> &dof_handler, 1104 IndexSet & dof_set) 1105 { 1106 // collect all the locally owned dofs 1107 dof_set = dof_handler.locally_owned_dofs(); 1108 1109 // add the DoF on the adjacent ghost cells to the IndexSet, cache them 1110 // in a set. need to check each dof manually because we can't be sure 1111 // that the dof range of locally_owned_dofs is really contiguous. 1112 std::vector<types::global_dof_index> dof_indices; 1113 std::set<types::global_dof_index> global_dof_indices; 1114 1115 typename DoFHandler<dim, spacedim>::active_cell_iterator 1116 cell = dof_handler.begin_active(), 1117 endc = dof_handler.end(); 1118 for (; cell != endc; ++cell) 1119 if (cell->is_locally_owned()) 1120 { 1121 dof_indices.resize(cell->get_fe().n_dofs_per_cell()); 1122 cell->get_dof_indices(dof_indices); 1123 1124 for (const types::global_dof_index dof_index : dof_indices) 1125 if (!dof_set.is_element(dof_index)) 1126 global_dof_indices.insert(dof_index); 1127 } 1128 1129 dof_set.add_indices(global_dof_indices.begin(), global_dof_indices.end()); 1130 1131 dof_set.compress(); 1132 } 1133 1134 1135 1136 template <int dim, int spacedim> 1137 void extract_locally_relevant_dofs(const DoFHandler<dim,spacedim> & dof_handler,IndexSet & dof_set)1138 extract_locally_relevant_dofs(const DoFHandler<dim, spacedim> &dof_handler, 1139 IndexSet & dof_set) 1140 { 1141 // collect all the locally owned dofs 1142 dof_set = dof_handler.locally_owned_dofs(); 1143 1144 // now add the DoF on the adjacent ghost cells to the IndexSet 1145 1146 // Note: For certain meshes (in particular in 3D and with many 1147 // processors), it is really necessary to cache intermediate data. After 1148 // trying several objects such as std::set, a vector that is always kept 1149 // sorted, and a vector that is initially unsorted and sorted once at the 1150 // end, the latter has been identified to provide the best performance. 1151 // Martin Kronbichler 1152 std::vector<types::global_dof_index> dof_indices; 1153 std::vector<types::global_dof_index> dofs_on_ghosts; 1154 1155 typename DoFHandler<dim, spacedim>::active_cell_iterator 1156 cell = dof_handler.begin_active(), 1157 endc = dof_handler.end(); 1158 for (; cell != endc; ++cell) 1159 if (cell->is_ghost()) 1160 { 1161 dof_indices.resize(cell->get_fe().n_dofs_per_cell()); 1162 cell->get_dof_indices(dof_indices); 1163 for (const auto dof_index : dof_indices) 1164 if (!dof_set.is_element(dof_index)) 1165 dofs_on_ghosts.push_back(dof_index); 1166 } 1167 1168 // sort, compress out duplicates, fill into index set 1169 std::sort(dofs_on_ghosts.begin(), dofs_on_ghosts.end()); 1170 dof_set.add_indices(dofs_on_ghosts.begin(), 1171 std::unique(dofs_on_ghosts.begin(), 1172 dofs_on_ghosts.end())); 1173 dof_set.compress(); 1174 } 1175 1176 1177 1178 template <int dim, int spacedim> 1179 void extract_locally_relevant_level_dofs(const DoFHandler<dim,spacedim> & dof_handler,const unsigned int level,IndexSet & dof_set)1180 extract_locally_relevant_level_dofs( 1181 const DoFHandler<dim, spacedim> &dof_handler, 1182 const unsigned int level, 1183 IndexSet & dof_set) 1184 { 1185 // collect all the locally owned dofs 1186 dof_set = dof_handler.locally_owned_mg_dofs(level); 1187 1188 // add the DoF on the adjacent ghost cells to the IndexSet 1189 1190 // Note: For certain meshes (in particular in 3D and with many 1191 // processors), it is really necessary to cache intermediate data. After 1192 // trying several objects such as std::set, a vector that is always kept 1193 // sorted, and a vector that is initially unsorted and sorted once at the 1194 // end, the latter has been identified to provide the best performance. 1195 // Martin Kronbichler 1196 std::vector<types::global_dof_index> dof_indices; 1197 std::vector<types::global_dof_index> dofs_on_ghosts; 1198 1199 typename DoFHandler<dim, spacedim>::cell_iterator cell = dof_handler.begin( 1200 level), 1201 endc = 1202 dof_handler.end(level); 1203 for (; cell != endc; ++cell) 1204 { 1205 const types::subdomain_id id = cell->level_subdomain_id(); 1206 1207 // skip artificial and own cells (only look at ghost cells) 1208 if (id == dof_handler.get_triangulation().locally_owned_subdomain() || 1209 id == numbers::artificial_subdomain_id) 1210 continue; 1211 1212 dof_indices.resize(cell->get_fe().n_dofs_per_cell()); 1213 cell->get_mg_dof_indices(dof_indices); 1214 for (const auto dof_index : dof_indices) 1215 if (!dof_set.is_element(dof_index)) 1216 dofs_on_ghosts.push_back(dof_index); 1217 } 1218 1219 // sort, compress out duplicates, fill into index set 1220 std::sort(dofs_on_ghosts.begin(), dofs_on_ghosts.end()); 1221 dof_set.add_indices(dofs_on_ghosts.begin(), 1222 std::unique(dofs_on_ghosts.begin(), 1223 dofs_on_ghosts.end())); 1224 1225 dof_set.compress(); 1226 } 1227 1228 1229 1230 template <int dim, int spacedim> 1231 void extract_constant_modes(const DoFHandler<dim,spacedim> & dof_handler,const ComponentMask & component_mask,std::vector<std::vector<bool>> & constant_modes)1232 extract_constant_modes(const DoFHandler<dim, spacedim> &dof_handler, 1233 const ComponentMask & component_mask, 1234 std::vector<std::vector<bool>> & constant_modes) 1235 { 1236 // If there are no locally owned DoFs, return with an empty 1237 // constant_modes object: 1238 if (dof_handler.n_locally_owned_dofs() == 0) 1239 { 1240 constant_modes = std::vector<std::vector<bool>>(0); 1241 return; 1242 } 1243 1244 const unsigned int n_components = dof_handler.get_fe(0).n_components(); 1245 Assert(component_mask.represents_n_components(n_components), 1246 ExcDimensionMismatch(n_components, component_mask.size())); 1247 1248 std::vector<unsigned char> dofs_by_component( 1249 dof_handler.n_locally_owned_dofs()); 1250 internal::get_component_association(dof_handler, 1251 component_mask, 1252 dofs_by_component); 1253 unsigned int n_selected_dofs = 0; 1254 for (unsigned int i = 0; i < n_components; ++i) 1255 if (component_mask[i] == true) 1256 n_selected_dofs += 1257 std::count(dofs_by_component.begin(), dofs_by_component.end(), i); 1258 1259 // Find local numbering within the selected components 1260 const IndexSet &locally_owned_dofs = dof_handler.locally_owned_dofs(); 1261 std::vector<unsigned int> component_numbering( 1262 locally_owned_dofs.n_elements(), numbers::invalid_unsigned_int); 1263 for (unsigned int i = 0, count = 0; i < locally_owned_dofs.n_elements(); 1264 ++i) 1265 if (component_mask[dofs_by_component[i]]) 1266 component_numbering[i] = count++; 1267 1268 // get the element constant modes and find a translation table between 1269 // index in the constant modes and the components. 1270 // 1271 // TODO: We might be able to extend this also for elements which do not 1272 // have the same constant modes, but that is messy... 1273 const dealii::hp::FECollection<dim, spacedim> &fe_collection = 1274 dof_handler.get_fe_collection(); 1275 std::vector<Table<2, bool>> element_constant_modes; 1276 std::vector<std::vector<std::pair<unsigned int, unsigned int>>> 1277 constant_mode_to_component_translation(n_components); 1278 unsigned int n_constant_modes = 0; 1279 for (unsigned int f = 0; f < fe_collection.size(); ++f) 1280 { 1281 std::pair<Table<2, bool>, std::vector<unsigned int>> data = 1282 fe_collection[f].get_constant_modes(); 1283 element_constant_modes.push_back(data.first); 1284 if (f == 0) 1285 for (unsigned int i = 0; i < data.second.size(); ++i) 1286 if (component_mask[data.second[i]]) 1287 constant_mode_to_component_translation[data.second[i]] 1288 .emplace_back(n_constant_modes++, i); 1289 AssertDimension(element_constant_modes.back().n_rows(), 1290 element_constant_modes[0].n_rows()); 1291 } 1292 1293 // First count the number of dofs in the current component. 1294 constant_modes.clear(); 1295 constant_modes.resize(n_constant_modes, 1296 std::vector<bool>(n_selected_dofs, false)); 1297 1298 // Loop over all owned cells and ask the element for the constant modes 1299 std::vector<types::global_dof_index> dof_indices; 1300 for (const auto &cell : dof_handler.active_cell_iterators()) 1301 if (cell->is_locally_owned()) 1302 { 1303 dof_indices.resize(cell->get_fe().n_dofs_per_cell()); 1304 cell->get_dof_indices(dof_indices); 1305 1306 for (unsigned int i = 0; i < dof_indices.size(); ++i) 1307 if (locally_owned_dofs.is_element(dof_indices[i])) 1308 { 1309 const unsigned int loc_index = 1310 locally_owned_dofs.index_within_set(dof_indices[i]); 1311 const unsigned int comp = dofs_by_component[loc_index]; 1312 if (component_mask[comp]) 1313 for (auto &indices : 1314 constant_mode_to_component_translation[comp]) 1315 constant_modes[indices 1316 .first][component_numbering[loc_index]] = 1317 element_constant_modes[cell->active_fe_index()]( 1318 indices.second, i); 1319 } 1320 } 1321 } 1322 1323 1324 1325 template <int dim, int spacedim> 1326 void get_active_fe_indices(const DoFHandler<dim,spacedim> & dof_handler,std::vector<unsigned int> & active_fe_indices)1327 get_active_fe_indices(const DoFHandler<dim, spacedim> &dof_handler, 1328 std::vector<unsigned int> & active_fe_indices) 1329 { 1330 AssertDimension(active_fe_indices.size(), 1331 dof_handler.get_triangulation().n_active_cells()); 1332 1333 typename DoFHandler<dim, spacedim>::active_cell_iterator 1334 cell = dof_handler.begin_active(), 1335 endc = dof_handler.end(); 1336 for (; cell != endc; ++cell) 1337 active_fe_indices[cell->active_cell_index()] = cell->active_fe_index(); 1338 } 1339 1340 template <int dim, int spacedim> 1341 std::vector<IndexSet> locally_owned_dofs_per_subdomain(const DoFHandler<dim,spacedim> & dof_handler)1342 locally_owned_dofs_per_subdomain(const DoFHandler<dim, spacedim> &dof_handler) 1343 { 1344 Assert(dof_handler.n_dofs() > 0, 1345 ExcMessage("The given DoFHandler has no DoFs.")); 1346 1347 // If the Triangulation is distributed, the only thing we can usefully 1348 // ask is for its locally owned subdomain 1349 Assert((dynamic_cast< 1350 const parallel::distributed::Triangulation<dim, spacedim> *>( 1351 &dof_handler.get_triangulation()) == nullptr), 1352 ExcMessage( 1353 "For parallel::distributed::Triangulation objects and " 1354 "associated DoF handler objects, asking for any information " 1355 "related to a subdomain other than the locally owned one does " 1356 "not make sense.")); 1357 1358 // The following is a random process (flip of a coin), thus should be called 1359 // once only. 1360 std::vector<dealii::types::subdomain_id> subdomain_association( 1361 dof_handler.n_dofs()); 1362 dealii::DoFTools::get_subdomain_association(dof_handler, 1363 subdomain_association); 1364 1365 // Figure out how many subdomain ids there are. 1366 // 1367 // if this is a parallel triangulation, then we can just ask the 1368 // triangulation for this. if this is a sequential triangulation, we loop 1369 // over all cells and take the largest subdomain_id value we find; the 1370 // number of subdomains is then the largest found value plus one. (we here 1371 // assume that all subdomain ids up to the largest are actually used; this 1372 // may not be true for a sequential triangulation where these values have 1373 // been set by hand and not in accordance with some MPI communicator; but 1374 // the function returns an array indexed starting at zero, so we need to 1375 // collect information for each subdomain index anyway, not just for the 1376 // used one.) 1377 const unsigned int n_subdomains = 1378 (dynamic_cast<const parallel::TriangulationBase<dim, spacedim> *>( 1379 &dof_handler.get_triangulation()) == nullptr ? 1380 [&dof_handler]() { 1381 unsigned int max_subdomain_id = 0; 1382 for (const auto &cell : dof_handler.active_cell_iterators()) 1383 max_subdomain_id = 1384 std::max(max_subdomain_id, cell->subdomain_id()); 1385 return max_subdomain_id + 1; 1386 }() : 1387 Utilities::MPI::n_mpi_processes( 1388 dynamic_cast<const parallel::TriangulationBase<dim, spacedim> *>( 1389 &dof_handler.get_triangulation()) 1390 ->get_communicator())); 1391 Assert(n_subdomains > *std::max_element(subdomain_association.begin(), 1392 subdomain_association.end()), 1393 ExcInternalError()); 1394 1395 std::vector<dealii::IndexSet> index_sets( 1396 n_subdomains, dealii::IndexSet(dof_handler.n_dofs())); 1397 1398 // loop over subdomain_association and populate IndexSet when a 1399 // change in subdomain ID is found 1400 dealii::types::global_dof_index i_min = 0; 1401 dealii::types::global_dof_index this_subdomain = subdomain_association[0]; 1402 1403 for (dealii::types::global_dof_index index = 1; 1404 index < subdomain_association.size(); 1405 ++index) 1406 { 1407 // found index different from the current one 1408 if (subdomain_association[index] != this_subdomain) 1409 { 1410 index_sets[this_subdomain].add_range(i_min, index); 1411 i_min = index; 1412 this_subdomain = subdomain_association[index]; 1413 } 1414 } 1415 1416 // the very last element is of different index 1417 if (i_min == subdomain_association.size() - 1) 1418 { 1419 index_sets[this_subdomain].add_index(i_min); 1420 } 1421 1422 // otherwise there are at least two different indices 1423 else 1424 { 1425 index_sets[this_subdomain].add_range(i_min, 1426 subdomain_association.size()); 1427 } 1428 1429 for (unsigned int i = 0; i < n_subdomains; i++) 1430 index_sets[i].compress(); 1431 1432 return index_sets; 1433 } 1434 1435 template <int dim, int spacedim> 1436 std::vector<IndexSet> locally_relevant_dofs_per_subdomain(const DoFHandler<dim,spacedim> & dof_handler)1437 locally_relevant_dofs_per_subdomain( 1438 const DoFHandler<dim, spacedim> &dof_handler) 1439 { 1440 // If the Triangulation is distributed, the only thing we can usefully 1441 // ask is for its locally owned subdomain 1442 Assert((dynamic_cast< 1443 const parallel::distributed::Triangulation<dim, spacedim> *>( 1444 &dof_handler.get_triangulation()) == nullptr), 1445 ExcMessage( 1446 "For parallel::distributed::Triangulation objects and " 1447 "associated DoF handler objects, asking for any information " 1448 "related to a subdomain other than the locally owned one does " 1449 "not make sense.")); 1450 1451 // Collect all the locally owned DoFs 1452 // Note: Even though the distribution of DoFs by the 1453 // locally_owned_dofs_per_subdomain function is pseudo-random, we will 1454 // collect all the DoFs on the subdomain and its layer cell. Therefore, the 1455 // random nature of this function does not play a role in the extraction of 1456 // the locally relevant DoFs 1457 std::vector<IndexSet> dof_set = 1458 locally_owned_dofs_per_subdomain(dof_handler); 1459 const dealii::types::subdomain_id n_subdomains = dof_set.size(); 1460 1461 // Add the DoFs on the adjacent (equivalent ghost) cells to the IndexSet, 1462 // cache them in a set. Need to check each DoF manually because we can't 1463 // be sure that the DoF range of locally_owned_dofs is really contiguous. 1464 for (dealii::types::subdomain_id subdomain_id = 0; 1465 subdomain_id < n_subdomains; 1466 ++subdomain_id) 1467 { 1468 // Extract the layer of cells around this subdomain 1469 std::function<bool( 1470 const typename DoFHandler<dim, spacedim>::active_cell_iterator &)> 1471 predicate = IteratorFilters::SubdomainEqualTo(subdomain_id); 1472 const std::vector< 1473 typename DoFHandler<dim, spacedim>::active_cell_iterator> 1474 active_halo_layer = 1475 GridTools::compute_active_cell_halo_layer(dof_handler, predicate); 1476 1477 // Extract DoFs associated with halo layer 1478 std::vector<types::global_dof_index> local_dof_indices; 1479 std::set<types::global_dof_index> subdomain_halo_global_dof_indices; 1480 for (typename std::vector< 1481 typename DoFHandler<dim, spacedim>::active_cell_iterator>:: 1482 const_iterator it_cell = active_halo_layer.begin(); 1483 it_cell != active_halo_layer.end(); 1484 ++it_cell) 1485 { 1486 const typename DoFHandler<dim, spacedim>::active_cell_iterator 1487 &cell = *it_cell; 1488 Assert( 1489 cell->subdomain_id() != subdomain_id, 1490 ExcMessage( 1491 "The subdomain ID of the halo cell should not match that of the vector entry.")); 1492 1493 local_dof_indices.resize(cell->get_fe().n_dofs_per_cell()); 1494 cell->get_dof_indices(local_dof_indices); 1495 1496 for (const types::global_dof_index local_dof_index : 1497 local_dof_indices) 1498 subdomain_halo_global_dof_indices.insert(local_dof_index); 1499 } 1500 1501 dof_set[subdomain_id].add_indices( 1502 subdomain_halo_global_dof_indices.begin(), 1503 subdomain_halo_global_dof_indices.end()); 1504 1505 dof_set[subdomain_id].compress(); 1506 } 1507 1508 return dof_set; 1509 } 1510 1511 template <int dim, int spacedim> 1512 void get_subdomain_association(const DoFHandler<dim,spacedim> & dof_handler,std::vector<types::subdomain_id> & subdomain_association)1513 get_subdomain_association( 1514 const DoFHandler<dim, spacedim> & dof_handler, 1515 std::vector<types::subdomain_id> &subdomain_association) 1516 { 1517 // if the Triangulation is distributed, the only thing we can usefully 1518 // ask is for its locally owned subdomain 1519 Assert((dynamic_cast< 1520 const parallel::distributed::Triangulation<dim, spacedim> *>( 1521 &dof_handler.get_triangulation()) == nullptr), 1522 ExcMessage( 1523 "For parallel::distributed::Triangulation objects and " 1524 "associated DoF handler objects, asking for any subdomain other " 1525 "than the locally owned one does not make sense.")); 1526 1527 Assert(subdomain_association.size() == dof_handler.n_dofs(), 1528 ExcDimensionMismatch(subdomain_association.size(), 1529 dof_handler.n_dofs())); 1530 1531 // catch an error that happened in some versions of the shared tria 1532 // distribute_dofs() function where we were trying to call this 1533 // function at a point in time when not all internal DoFHandler 1534 // structures were quite set up yet. 1535 Assert(dof_handler.n_dofs() > 0, ExcInternalError()); 1536 1537 // In case this function is executed with parallel::shared::Triangulation 1538 // with possibly artificial cells, we need to take "true" subdomain IDs 1539 // (i.e. without artificial cells). Otherwise we are good to use 1540 // subdomain_id as stored in cell->subdomain_id(). 1541 std::vector<types::subdomain_id> cell_owners( 1542 dof_handler.get_triangulation().n_active_cells()); 1543 if (const parallel::shared::Triangulation<dim, spacedim> *tr = 1544 (dynamic_cast<const parallel::shared::Triangulation<dim, spacedim> *>( 1545 &dof_handler.get_triangulation()))) 1546 { 1547 cell_owners = tr->get_true_subdomain_ids_of_cells(); 1548 Assert(tr->get_true_subdomain_ids_of_cells().size() == 1549 tr->n_active_cells(), 1550 ExcInternalError()); 1551 } 1552 else 1553 { 1554 for (typename DoFHandler<dim, spacedim>::active_cell_iterator cell = 1555 dof_handler.begin_active(); 1556 cell != dof_handler.end(); 1557 cell++) 1558 if (cell->is_locally_owned()) 1559 cell_owners[cell->active_cell_index()] = cell->subdomain_id(); 1560 } 1561 1562 // preset all values by an invalid value 1563 std::fill_n(subdomain_association.begin(), 1564 dof_handler.n_dofs(), 1565 numbers::invalid_subdomain_id); 1566 1567 std::vector<types::global_dof_index> local_dof_indices; 1568 local_dof_indices.reserve( 1569 dof_handler.get_fe_collection().max_dofs_per_cell()); 1570 1571 // loop over all cells and record which subdomain a DoF belongs to. 1572 // give to the smaller subdomain_id in case it is on an interface 1573 typename DoFHandler<dim, spacedim>::active_cell_iterator 1574 cell = dof_handler.begin_active(), 1575 endc = dof_handler.end(); 1576 for (; cell != endc; ++cell) 1577 { 1578 const types::subdomain_id subdomain_id = 1579 cell_owners[cell->active_cell_index()]; 1580 const unsigned int dofs_per_cell = cell->get_fe().n_dofs_per_cell(); 1581 local_dof_indices.resize(dofs_per_cell); 1582 cell->get_dof_indices(local_dof_indices); 1583 1584 // set subdomain ids. if dofs already have their values set then 1585 // they must be on partition interfaces. in that case assign them 1586 // to either the previous association or the current processor 1587 // with the smaller subdomain id. 1588 for (unsigned int i = 0; i < dofs_per_cell; ++i) 1589 if (subdomain_association[local_dof_indices[i]] == 1590 numbers::invalid_subdomain_id) 1591 subdomain_association[local_dof_indices[i]] = subdomain_id; 1592 else if (subdomain_association[local_dof_indices[i]] > subdomain_id) 1593 { 1594 subdomain_association[local_dof_indices[i]] = subdomain_id; 1595 } 1596 } 1597 1598 Assert(std::find(subdomain_association.begin(), 1599 subdomain_association.end(), 1600 numbers::invalid_subdomain_id) == 1601 subdomain_association.end(), 1602 ExcInternalError()); 1603 } 1604 1605 1606 1607 template <int dim, int spacedim> 1608 unsigned int count_dofs_with_subdomain_association(const DoFHandler<dim,spacedim> & dof_handler,const types::subdomain_id subdomain)1609 count_dofs_with_subdomain_association( 1610 const DoFHandler<dim, spacedim> &dof_handler, 1611 const types::subdomain_id subdomain) 1612 { 1613 std::vector<types::subdomain_id> subdomain_association( 1614 dof_handler.n_dofs()); 1615 get_subdomain_association(dof_handler, subdomain_association); 1616 1617 return std::count(subdomain_association.begin(), 1618 subdomain_association.end(), 1619 subdomain); 1620 } 1621 1622 1623 1624 template <int dim, int spacedim> 1625 IndexSet dof_indices_with_subdomain_association(const DoFHandler<dim,spacedim> & dof_handler,const types::subdomain_id subdomain)1626 dof_indices_with_subdomain_association( 1627 const DoFHandler<dim, spacedim> &dof_handler, 1628 const types::subdomain_id subdomain) 1629 { 1630 // If we have a distributed::Triangulation only allow locally_owned 1631 // subdomain. 1632 Assert((dof_handler.get_triangulation().locally_owned_subdomain() == 1633 numbers::invalid_subdomain_id) || 1634 (subdomain == 1635 dof_handler.get_triangulation().locally_owned_subdomain()), 1636 ExcMessage( 1637 "For parallel::distributed::Triangulation objects and " 1638 "associated DoF handler objects, asking for any subdomain other " 1639 "than the locally owned one does not make sense.")); 1640 1641 IndexSet index_set(dof_handler.n_dofs()); 1642 1643 std::vector<types::global_dof_index> local_dof_indices; 1644 local_dof_indices.reserve( 1645 dof_handler.get_fe_collection().max_dofs_per_cell()); 1646 1647 // first generate an unsorted list of all indices which we fill from 1648 // the back. could also insert them directly into the IndexSet, but 1649 // that inserts indices in the middle, which is an O(n^2) algorithm and 1650 // hence too expensive. Could also use std::set, but that is in general 1651 // more expensive than a vector 1652 std::vector<types::global_dof_index> subdomain_indices; 1653 1654 typename DoFHandler<dim, spacedim>::active_cell_iterator 1655 cell = dof_handler.begin_active(), 1656 endc = dof_handler.end(); 1657 for (; cell != endc; ++cell) 1658 if ((cell->is_artificial() == false) && 1659 (cell->subdomain_id() == subdomain)) 1660 { 1661 const unsigned int dofs_per_cell = cell->get_fe().n_dofs_per_cell(); 1662 local_dof_indices.resize(dofs_per_cell); 1663 cell->get_dof_indices(local_dof_indices); 1664 subdomain_indices.insert(subdomain_indices.end(), 1665 local_dof_indices.begin(), 1666 local_dof_indices.end()); 1667 } 1668 // sort indices and remove duplicates 1669 std::sort(subdomain_indices.begin(), subdomain_indices.end()); 1670 subdomain_indices.erase(std::unique(subdomain_indices.begin(), 1671 subdomain_indices.end()), 1672 subdomain_indices.end()); 1673 1674 // insert into IndexSet 1675 index_set.add_indices(subdomain_indices.begin(), subdomain_indices.end()); 1676 index_set.compress(); 1677 1678 return index_set; 1679 } 1680 1681 1682 1683 template <int dim, int spacedim> 1684 void count_dofs_with_subdomain_association(const DoFHandler<dim,spacedim> & dof_handler,const types::subdomain_id subdomain,std::vector<unsigned int> & n_dofs_on_subdomain)1685 count_dofs_with_subdomain_association( 1686 const DoFHandler<dim, spacedim> &dof_handler, 1687 const types::subdomain_id subdomain, 1688 std::vector<unsigned int> & n_dofs_on_subdomain) 1689 { 1690 Assert(n_dofs_on_subdomain.size() == dof_handler.get_fe(0).n_components(), 1691 ExcDimensionMismatch(n_dofs_on_subdomain.size(), 1692 dof_handler.get_fe(0).n_components())); 1693 std::fill(n_dofs_on_subdomain.begin(), n_dofs_on_subdomain.end(), 0); 1694 1695 // Make sure there are at least some cells with this subdomain id 1696 Assert(std::any_of( 1697 dof_handler.begin_active(), 1698 typename DoFHandler<dim, spacedim>::active_cell_iterator{ 1699 dof_handler.end()}, 1700 [subdomain]( 1701 const typename DoFHandler<dim, spacedim>::cell_accessor &cell) { 1702 return cell.subdomain_id() == subdomain; 1703 }), 1704 ExcMessage("There are no cells for the given subdomain!")); 1705 1706 std::vector<types::subdomain_id> subdomain_association( 1707 dof_handler.n_dofs()); 1708 get_subdomain_association(dof_handler, subdomain_association); 1709 1710 std::vector<unsigned char> component_association(dof_handler.n_dofs()); 1711 internal::get_component_association(dof_handler, 1712 std::vector<bool>(), 1713 component_association); 1714 1715 for (unsigned int c = 0; c < dof_handler.get_fe(0).n_components(); ++c) 1716 { 1717 for (types::global_dof_index i = 0; i < dof_handler.n_dofs(); ++i) 1718 if ((subdomain_association[i] == subdomain) && 1719 (component_association[i] == static_cast<unsigned char>(c))) 1720 ++n_dofs_on_subdomain[c]; 1721 } 1722 } 1723 1724 1725 1726 namespace internal 1727 { 1728 // TODO: why is this function so complicated? It would be nice to have 1729 // comments that explain why we can't just loop over all components and 1730 // count the entries in dofs_by_component that have this component's 1731 // index 1732 template <int dim, int spacedim> 1733 void resolve_components(const FiniteElement<dim,spacedim> & fe,const std::vector<unsigned char> & dofs_by_component,const std::vector<unsigned int> & target_component,const bool only_once,std::vector<types::global_dof_index> & dofs_per_component,unsigned int & component)1734 resolve_components(const FiniteElement<dim, spacedim> & fe, 1735 const std::vector<unsigned char> & dofs_by_component, 1736 const std::vector<unsigned int> & target_component, 1737 const bool only_once, 1738 std::vector<types::global_dof_index> &dofs_per_component, 1739 unsigned int & component) 1740 { 1741 for (unsigned int b = 0; b < fe.n_base_elements(); ++b) 1742 { 1743 const FiniteElement<dim, spacedim> &base = fe.base_element(b); 1744 // Dimension of base element 1745 unsigned int d = base.n_components(); 1746 1747 for (unsigned int m = 0; m < fe.element_multiplicity(b); ++m) 1748 { 1749 if (base.n_base_elements() > 1) 1750 resolve_components(base, 1751 dofs_by_component, 1752 target_component, 1753 only_once, 1754 dofs_per_component, 1755 component); 1756 else 1757 { 1758 for (unsigned int dd = 0; dd < d; ++dd, ++component) 1759 dofs_per_component[target_component[component]] += 1760 std::count(dofs_by_component.begin(), 1761 dofs_by_component.end(), 1762 component); 1763 1764 // if we have non-primitive FEs and want all components 1765 // to show the number of dofs, need to copy the result to 1766 // those components 1767 if (!base.is_primitive() && !only_once) 1768 for (unsigned int dd = 1; dd < d; ++dd) 1769 dofs_per_component[target_component[component - d + dd]] = 1770 dofs_per_component[target_component[component - d]]; 1771 } 1772 } 1773 } 1774 } 1775 1776 1777 template <int dim, int spacedim> 1778 void resolve_components(const hp::FECollection<dim,spacedim> & fe_collection,const std::vector<unsigned char> & dofs_by_component,const std::vector<unsigned int> & target_component,const bool only_once,std::vector<types::global_dof_index> & dofs_per_component,unsigned int & component)1779 resolve_components(const hp::FECollection<dim, spacedim> &fe_collection, 1780 const std::vector<unsigned char> & dofs_by_component, 1781 const std::vector<unsigned int> & target_component, 1782 const bool only_once, 1783 std::vector<types::global_dof_index> &dofs_per_component, 1784 unsigned int & component) 1785 { 1786 // assert that all elements in the collection have the same structure 1787 // (base elements and multiplicity, components per base element) and 1788 // then simply call the function above 1789 for (unsigned int fe = 1; fe < fe_collection.size(); ++fe) 1790 { 1791 Assert(fe_collection[fe].n_components() == 1792 fe_collection[0].n_components(), 1793 ExcNotImplemented()); 1794 Assert(fe_collection[fe].n_base_elements() == 1795 fe_collection[0].n_base_elements(), 1796 ExcNotImplemented()); 1797 for (unsigned int b = 0; b < fe_collection[0].n_base_elements(); ++b) 1798 { 1799 Assert(fe_collection[fe].base_element(b).n_components() == 1800 fe_collection[0].base_element(b).n_components(), 1801 ExcNotImplemented()); 1802 Assert(fe_collection[fe].base_element(b).n_base_elements() == 1803 fe_collection[0].base_element(b).n_base_elements(), 1804 ExcNotImplemented()); 1805 } 1806 } 1807 1808 resolve_components(fe_collection[0], 1809 dofs_by_component, 1810 target_component, 1811 only_once, 1812 dofs_per_component, 1813 component); 1814 } 1815 } // namespace internal 1816 1817 1818 1819 namespace internal 1820 { 1821 namespace 1822 { 1823 /** 1824 * Return true if the given element is primitive. 1825 */ 1826 template <int dim, int spacedim> 1827 bool all_elements_are_primitive(const FiniteElement<dim,spacedim> & fe)1828 all_elements_are_primitive(const FiniteElement<dim, spacedim> &fe) 1829 { 1830 return fe.is_primitive(); 1831 } 1832 1833 1834 /** 1835 * Return true if each element of the given element collection is 1836 * primitive. 1837 */ 1838 template <int dim, int spacedim> 1839 bool all_elements_are_primitive(const dealii::hp::FECollection<dim,spacedim> & fe_collection)1840 all_elements_are_primitive( 1841 const dealii::hp::FECollection<dim, spacedim> &fe_collection) 1842 { 1843 for (unsigned int i = 0; i < fe_collection.size(); ++i) 1844 if (fe_collection[i].is_primitive() == false) 1845 return false; 1846 1847 return true; 1848 } 1849 } // namespace 1850 } // namespace internal 1851 1852 1853 1854 // deprecated function 1855 template <int dim, int spacedim> 1856 void count_dofs_per_component(const DoFHandler<dim,spacedim> & dof_handler,std::vector<types::global_dof_index> & dofs_per_component,const bool only_once,const std::vector<unsigned int> & target_component)1857 count_dofs_per_component( 1858 const DoFHandler<dim, spacedim> & dof_handler, 1859 std::vector<types::global_dof_index> &dofs_per_component, 1860 const bool only_once, 1861 const std::vector<unsigned int> & target_component) 1862 { 1863 dofs_per_component = 1864 count_dofs_per_fe_component(dof_handler, only_once, target_component); 1865 } 1866 1867 1868 1869 template <int dim, int spacedim> 1870 std::vector<types::global_dof_index> count_dofs_per_fe_component(const DoFHandler<dim,spacedim> & dof_handler,const bool only_once,const std::vector<unsigned int> & target_component_)1871 count_dofs_per_fe_component( 1872 const DoFHandler<dim, spacedim> &dof_handler, 1873 const bool only_once, 1874 const std::vector<unsigned int> &target_component_) 1875 { 1876 const unsigned int n_components = dof_handler.get_fe(0).n_components(); 1877 1878 // If the empty vector was given as default argument, set up this 1879 // vector as identity. 1880 std::vector<unsigned int> target_component = target_component_; 1881 if (target_component.size() == 0) 1882 { 1883 target_component.resize(n_components); 1884 for (unsigned int i = 0; i < n_components; ++i) 1885 target_component[i] = i; 1886 } 1887 else 1888 Assert(target_component.size() == n_components, 1889 ExcDimensionMismatch(target_component.size(), n_components)); 1890 1891 1892 const unsigned int max_component = 1893 *std::max_element(target_component.begin(), target_component.end()); 1894 const unsigned int n_target_components = max_component + 1; 1895 1896 std::vector<types::global_dof_index> dofs_per_component( 1897 n_target_components, types::global_dof_index(0)); 1898 1899 // special case for only one component. treat this first since it does 1900 // not require any computations 1901 if (n_components == 1) 1902 { 1903 dofs_per_component[0] = dof_handler.n_locally_owned_dofs(); 1904 return dofs_per_component; 1905 } 1906 1907 1908 // otherwise determine the number of dofs in each component separately. 1909 // do so in parallel 1910 std::vector<unsigned char> dofs_by_component( 1911 dof_handler.n_locally_owned_dofs()); 1912 internal::get_component_association(dof_handler, 1913 ComponentMask(), 1914 dofs_by_component); 1915 1916 // next count what we got 1917 unsigned int component = 0; 1918 internal::resolve_components(dof_handler.get_fe_collection(), 1919 dofs_by_component, 1920 target_component, 1921 only_once, 1922 dofs_per_component, 1923 component); 1924 Assert(n_components == component, ExcInternalError()); 1925 1926 // finally sanity check. this is only valid if the finite element is 1927 // actually primitive, so exclude other elements from this 1928 Assert((internal::all_elements_are_primitive( 1929 dof_handler.get_fe_collection()) == false) || 1930 (std::accumulate(dofs_per_component.begin(), 1931 dofs_per_component.end(), 1932 types::global_dof_index(0)) == 1933 dof_handler.n_locally_owned_dofs()), 1934 ExcInternalError()); 1935 1936 // reduce information from all CPUs 1937 #ifdef DEAL_II_WITH_MPI 1938 1939 if (const parallel::TriangulationBase<dim, spacedim> *tria = 1940 (dynamic_cast<const parallel::TriangulationBase<dim, spacedim> *>( 1941 &dof_handler.get_triangulation()))) 1942 { 1943 std::vector<types::global_dof_index> local_dof_count = 1944 dofs_per_component; 1945 1946 const int ierr = MPI_Allreduce(local_dof_count.data(), 1947 dofs_per_component.data(), 1948 n_target_components, 1949 DEAL_II_DOF_INDEX_MPI_TYPE, 1950 MPI_SUM, 1951 tria->get_communicator()); 1952 AssertThrowMPI(ierr); 1953 } 1954 #endif 1955 1956 return dofs_per_component; 1957 } 1958 1959 1960 1961 // deprecated function 1962 template <int dim, int spacedim> 1963 void count_dofs_per_block(const DoFHandler<dim,spacedim> & dof_handler,std::vector<types::global_dof_index> & dofs_per_block,const std::vector<unsigned int> & target_block)1964 count_dofs_per_block(const DoFHandler<dim, spacedim> & dof_handler, 1965 std::vector<types::global_dof_index> &dofs_per_block, 1966 const std::vector<unsigned int> & target_block) 1967 { 1968 dofs_per_block = count_dofs_per_fe_block(dof_handler, target_block); 1969 } 1970 1971 1972 1973 template <int dim, int spacedim> 1974 std::vector<types::global_dof_index> count_dofs_per_fe_block(const DoFHandler<dim,spacedim> & dof_handler,const std::vector<unsigned int> & target_block_)1975 count_dofs_per_fe_block(const DoFHandler<dim, spacedim> &dof_handler, 1976 const std::vector<unsigned int> &target_block_) 1977 { 1978 const dealii::hp::FECollection<dim, spacedim> &fe_collection = 1979 dof_handler.get_fe_collection(); 1980 Assert(fe_collection.size() < 256, ExcNotImplemented()); 1981 1982 // If the empty vector for target_block(e.g., as default argument), then 1983 // set up this vector as identity. We do this set up with the first 1984 // element of the collection, but the whole thing can only work if 1985 // all elements have the same number of blocks anyway -- so check 1986 // that right after 1987 const unsigned int n_blocks = fe_collection[0].n_blocks(); 1988 1989 std::vector<unsigned int> target_block = target_block_; 1990 if (target_block.size() == 0) 1991 { 1992 target_block.resize(fe_collection[0].n_blocks()); 1993 for (unsigned int i = 0; i < n_blocks; ++i) 1994 target_block[i] = i; 1995 } 1996 else 1997 Assert(target_block.size() == n_blocks, 1998 ExcDimensionMismatch(target_block.size(), n_blocks)); 1999 for (unsigned int f = 1; f < fe_collection.size(); ++f) 2000 Assert(fe_collection[0].n_blocks() == fe_collection[f].n_blocks(), 2001 ExcMessage("This function can only work if all elements in a " 2002 "collection have the same number of blocks.")); 2003 2004 // special case for only one block. treat this first since it does 2005 // not require any computations 2006 if (n_blocks == 1) 2007 { 2008 std::vector<types::global_dof_index> dofs_per_block(1); 2009 dofs_per_block[0] = dof_handler.n_dofs(); 2010 return dofs_per_block; 2011 } 2012 2013 // Otherwise set up the right-sized object and start working 2014 const unsigned int max_block = 2015 *std::max_element(target_block.begin(), target_block.end()); 2016 const unsigned int n_target_blocks = max_block + 1; 2017 2018 std::vector<types::global_dof_index> dofs_per_block(n_target_blocks); 2019 2020 // Loop over the elements of the collection, but really only consider 2021 // the last element (see #9271) 2022 for (unsigned int this_fe = fe_collection.size() - 1; 2023 this_fe < fe_collection.size(); 2024 ++this_fe) 2025 { 2026 const FiniteElement<dim, spacedim> &fe = fe_collection[this_fe]; 2027 2028 std::vector<unsigned char> dofs_by_block( 2029 dof_handler.n_locally_owned_dofs()); 2030 internal::get_block_association(dof_handler, dofs_by_block); 2031 2032 // next count what we got 2033 for (unsigned int block = 0; block < fe.n_blocks(); ++block) 2034 dofs_per_block[target_block[block]] += 2035 std::count(dofs_by_block.begin(), dofs_by_block.end(), block); 2036 2037 #ifdef DEAL_II_WITH_MPI 2038 // if we are working on a parallel mesh, we now need to collect 2039 // this information from all processors 2040 if (const parallel::TriangulationBase<dim, spacedim> *tria = 2041 (dynamic_cast<const parallel::TriangulationBase<dim, spacedim> *>( 2042 &dof_handler.get_triangulation()))) 2043 { 2044 std::vector<types::global_dof_index> local_dof_count = 2045 dofs_per_block; 2046 const int ierr = MPI_Allreduce(local_dof_count.data(), 2047 dofs_per_block.data(), 2048 n_target_blocks, 2049 DEAL_II_DOF_INDEX_MPI_TYPE, 2050 MPI_SUM, 2051 tria->get_communicator()); 2052 AssertThrowMPI(ierr); 2053 } 2054 #endif 2055 } 2056 2057 return dofs_per_block; 2058 } 2059 2060 2061 2062 template <int dim, int spacedim> 2063 void map_dof_to_boundary_indices(const DoFHandler<dim,spacedim> & dof_handler,std::vector<types::global_dof_index> & mapping)2064 map_dof_to_boundary_indices(const DoFHandler<dim, spacedim> & dof_handler, 2065 std::vector<types::global_dof_index> &mapping) 2066 { 2067 mapping.clear(); 2068 mapping.insert(mapping.end(), 2069 dof_handler.n_dofs(), 2070 numbers::invalid_dof_index); 2071 2072 std::vector<types::global_dof_index> dofs_on_face; 2073 dofs_on_face.reserve(dof_handler.get_fe_collection().max_dofs_per_face()); 2074 types::global_dof_index next_boundary_index = 0; 2075 2076 // now loop over all cells and check whether their faces are at the 2077 // boundary. note that we need not take special care of single lines 2078 // being at the boundary (using @p{cell->has_boundary_lines}), since we 2079 // do not support boundaries of dimension dim-2, and so every isolated 2080 // boundary line is also part of a boundary face which we will be 2081 // visiting sooner or later 2082 typename DoFHandler<dim, spacedim>::active_cell_iterator 2083 cell = dof_handler.begin_active(), 2084 endc = dof_handler.end(); 2085 for (; cell != endc; ++cell) 2086 for (const unsigned int f : cell->face_indices()) 2087 if (cell->at_boundary(f)) 2088 { 2089 const unsigned int dofs_per_face = cell->get_fe().n_dofs_per_face(); 2090 dofs_on_face.resize(dofs_per_face); 2091 cell->face(f)->get_dof_indices(dofs_on_face, 2092 cell->active_fe_index()); 2093 for (unsigned int i = 0; i < dofs_per_face; ++i) 2094 if (mapping[dofs_on_face[i]] == numbers::invalid_dof_index) 2095 mapping[dofs_on_face[i]] = next_boundary_index++; 2096 } 2097 2098 AssertDimension(next_boundary_index, dof_handler.n_boundary_dofs()); 2099 } 2100 2101 2102 2103 template <int dim, int spacedim> 2104 void map_dof_to_boundary_indices(const DoFHandler<dim,spacedim> & dof_handler,const std::set<types::boundary_id> & boundary_ids,std::vector<types::global_dof_index> & mapping)2105 map_dof_to_boundary_indices(const DoFHandler<dim, spacedim> & dof_handler, 2106 const std::set<types::boundary_id> &boundary_ids, 2107 std::vector<types::global_dof_index> &mapping) 2108 { 2109 Assert(boundary_ids.find(numbers::internal_face_boundary_id) == 2110 boundary_ids.end(), 2111 ExcInvalidBoundaryIndicator()); 2112 2113 mapping.clear(); 2114 mapping.insert(mapping.end(), 2115 dof_handler.n_dofs(), 2116 numbers::invalid_dof_index); 2117 2118 // return if there is nothing to do 2119 if (boundary_ids.size() == 0) 2120 return; 2121 2122 std::vector<types::global_dof_index> dofs_on_face; 2123 dofs_on_face.reserve(dof_handler.get_fe_collection().max_dofs_per_face()); 2124 types::global_dof_index next_boundary_index = 0; 2125 2126 typename DoFHandler<dim, spacedim>::active_cell_iterator 2127 cell = dof_handler.begin_active(), 2128 endc = dof_handler.end(); 2129 for (; cell != endc; ++cell) 2130 for (const unsigned int f : cell->face_indices()) 2131 if (boundary_ids.find(cell->face(f)->boundary_id()) != 2132 boundary_ids.end()) 2133 { 2134 const unsigned int dofs_per_face = cell->get_fe().n_dofs_per_face(); 2135 dofs_on_face.resize(dofs_per_face); 2136 cell->face(f)->get_dof_indices(dofs_on_face, 2137 cell->active_fe_index()); 2138 for (unsigned int i = 0; i < dofs_per_face; ++i) 2139 if (mapping[dofs_on_face[i]] == numbers::invalid_dof_index) 2140 mapping[dofs_on_face[i]] = next_boundary_index++; 2141 } 2142 2143 AssertDimension(next_boundary_index, 2144 dof_handler.n_boundary_dofs(boundary_ids)); 2145 } 2146 2147 namespace internal 2148 { 2149 namespace 2150 { 2151 template <int dim, int spacedim> 2152 void map_dofs_to_support_points(const hp::MappingCollection<dim,spacedim> & mapping,const DoFHandler<dim,spacedim> & dof_handler,std::map<types::global_dof_index,Point<spacedim>> & support_points,const ComponentMask & in_mask)2153 map_dofs_to_support_points( 2154 const hp::MappingCollection<dim, spacedim> & mapping, 2155 const DoFHandler<dim, spacedim> & dof_handler, 2156 std::map<types::global_dof_index, Point<spacedim>> &support_points, 2157 const ComponentMask & in_mask) 2158 { 2159 const hp::FECollection<dim, spacedim> &fe_collection = 2160 dof_handler.get_fe_collection(); 2161 hp::QCollection<dim> q_coll_dummy; 2162 2163 for (unsigned int fe_index = 0; fe_index < fe_collection.size(); 2164 ++fe_index) 2165 { 2166 // check whether every fe in the collection has support points 2167 Assert(fe_collection[fe_index].has_support_points(), 2168 typename FiniteElement<dim>::ExcFEHasNoSupportPoints()); 2169 q_coll_dummy.push_back(Quadrature<dim>( 2170 fe_collection[fe_index].get_unit_support_points())); 2171 } 2172 2173 // Take care of components 2174 const ComponentMask mask = 2175 (in_mask.size() == 0 ? 2176 ComponentMask(fe_collection.n_components(), true) : 2177 in_mask); 2178 2179 // Now loop over all cells and enquire the support points on each 2180 // of these. we use dummy quadrature formulas where the quadrature 2181 // points are located at the unit support points to enquire the 2182 // location of the support points in real space. 2183 // 2184 // The weights of the quadrature rule have been set to invalid 2185 // values by the used constructor. 2186 hp::FEValues<dim, spacedim> hp_fe_values(mapping, 2187 fe_collection, 2188 q_coll_dummy, 2189 update_quadrature_points); 2190 typename DoFHandler<dim, spacedim>::active_cell_iterator 2191 cell = dof_handler.begin_active(), 2192 endc = dof_handler.end(); 2193 2194 std::vector<types::global_dof_index> local_dof_indices; 2195 for (; cell != endc; ++cell) 2196 // only work on locally relevant cells 2197 if (cell->is_artificial() == false) 2198 { 2199 hp_fe_values.reinit(cell); 2200 const FEValues<dim, spacedim> &fe_values = 2201 hp_fe_values.get_present_fe_values(); 2202 2203 local_dof_indices.resize(cell->get_fe().n_dofs_per_cell()); 2204 cell->get_dof_indices(local_dof_indices); 2205 2206 const std::vector<Point<spacedim>> &points = 2207 fe_values.get_quadrature_points(); 2208 for (unsigned int i = 0; i < cell->get_fe().n_dofs_per_cell(); 2209 ++i) 2210 { 2211 const unsigned int dof_comp = 2212 cell->get_fe().system_to_component_index(i).first; 2213 2214 // insert the values into the map if it is a valid component 2215 if (mask[dof_comp]) 2216 support_points[local_dof_indices[i]] = points[i]; 2217 } 2218 } 2219 } 2220 2221 2222 template <int dim, int spacedim> 2223 void map_dofs_to_support_points(const hp::MappingCollection<dim,spacedim> & mapping,const DoFHandler<dim,spacedim> & dof_handler,std::vector<Point<spacedim>> & support_points,const ComponentMask & mask)2224 map_dofs_to_support_points( 2225 const hp::MappingCollection<dim, spacedim> &mapping, 2226 const DoFHandler<dim, spacedim> & dof_handler, 2227 std::vector<Point<spacedim>> & support_points, 2228 const ComponentMask & mask) 2229 { 2230 // get the data in the form of the map as above 2231 std::map<types::global_dof_index, Point<spacedim>> x_support_points; 2232 map_dofs_to_support_points(mapping, 2233 dof_handler, 2234 x_support_points, 2235 mask); 2236 2237 // now convert from the map to the linear vector. make sure every 2238 // entry really appeared in the map 2239 for (types::global_dof_index i = 0; i < dof_handler.n_dofs(); ++i) 2240 { 2241 Assert(x_support_points.find(i) != x_support_points.end(), 2242 ExcInternalError()); 2243 2244 support_points[i] = x_support_points[i]; 2245 } 2246 } 2247 } // namespace 2248 } // namespace internal 2249 2250 template <int dim, int spacedim> 2251 void map_dofs_to_support_points(const Mapping<dim,spacedim> & mapping,const DoFHandler<dim,spacedim> & dof_handler,std::vector<Point<spacedim>> & support_points,const ComponentMask & mask)2252 map_dofs_to_support_points(const Mapping<dim, spacedim> & mapping, 2253 const DoFHandler<dim, spacedim> &dof_handler, 2254 std::vector<Point<spacedim>> & support_points, 2255 const ComponentMask & mask) 2256 { 2257 AssertDimension(support_points.size(), dof_handler.n_dofs()); 2258 Assert((dynamic_cast< 2259 const parallel::distributed::Triangulation<dim, spacedim> *>( 2260 &dof_handler.get_triangulation()) == nullptr), 2261 ExcMessage( 2262 "This function can not be used with distributed triangulations. " 2263 "See the documentation for more information.")); 2264 2265 // Let the internal function do all the work, just make sure that it 2266 // gets a MappingCollection 2267 const hp::MappingCollection<dim, spacedim> mapping_collection(mapping); 2268 2269 internal::map_dofs_to_support_points(mapping_collection, 2270 dof_handler, 2271 support_points, 2272 mask); 2273 } 2274 2275 2276 template <int dim, int spacedim> 2277 void map_dofs_to_support_points(const hp::MappingCollection<dim,spacedim> & mapping,const DoFHandler<dim,spacedim> & dof_handler,std::vector<Point<spacedim>> & support_points,const ComponentMask & mask)2278 map_dofs_to_support_points( 2279 const hp::MappingCollection<dim, spacedim> &mapping, 2280 const DoFHandler<dim, spacedim> & dof_handler, 2281 std::vector<Point<spacedim>> & support_points, 2282 const ComponentMask & mask) 2283 { 2284 AssertDimension(support_points.size(), dof_handler.n_dofs()); 2285 Assert((dynamic_cast< 2286 const parallel::distributed::Triangulation<dim, spacedim> *>( 2287 &dof_handler.get_triangulation()) == nullptr), 2288 ExcMessage( 2289 "This function can not be used with distributed triangulations. " 2290 "See the documentation for more information.")); 2291 2292 // Let the internal function do all the work, just make sure that it 2293 // gets a MappingCollection 2294 internal::map_dofs_to_support_points(mapping, 2295 dof_handler, 2296 support_points, 2297 mask); 2298 } 2299 2300 2301 template <int dim, int spacedim> 2302 void map_dofs_to_support_points(const Mapping<dim,spacedim> & mapping,const DoFHandler<dim,spacedim> & dof_handler,std::map<types::global_dof_index,Point<spacedim>> & support_points,const ComponentMask & mask)2303 map_dofs_to_support_points( 2304 const Mapping<dim, spacedim> & mapping, 2305 const DoFHandler<dim, spacedim> & dof_handler, 2306 std::map<types::global_dof_index, Point<spacedim>> &support_points, 2307 const ComponentMask & mask) 2308 { 2309 support_points.clear(); 2310 2311 // Let the internal function do all the work, just make sure that it 2312 // gets a MappingCollection 2313 const hp::MappingCollection<dim, spacedim> mapping_collection(mapping); 2314 2315 internal::map_dofs_to_support_points(mapping_collection, 2316 dof_handler, 2317 support_points, 2318 mask); 2319 } 2320 2321 2322 template <int dim, int spacedim> 2323 void map_dofs_to_support_points(const hp::MappingCollection<dim,spacedim> & mapping,const DoFHandler<dim,spacedim> & dof_handler,std::map<types::global_dof_index,Point<spacedim>> & support_points,const ComponentMask & mask)2324 map_dofs_to_support_points( 2325 const hp::MappingCollection<dim, spacedim> & mapping, 2326 const DoFHandler<dim, spacedim> & dof_handler, 2327 std::map<types::global_dof_index, Point<spacedim>> &support_points, 2328 const ComponentMask & mask) 2329 { 2330 support_points.clear(); 2331 2332 // Let the internal function do all the work, just make sure that it 2333 // gets a MappingCollection 2334 internal::map_dofs_to_support_points(mapping, 2335 dof_handler, 2336 support_points, 2337 mask); 2338 } 2339 2340 template <int spacedim> 2341 void write_gnuplot_dof_support_point_info(std::ostream & out,const std::map<types::global_dof_index,Point<spacedim>> & support_points)2342 write_gnuplot_dof_support_point_info( 2343 std::ostream & out, 2344 const std::map<types::global_dof_index, Point<spacedim>> &support_points) 2345 { 2346 AssertThrow(out, ExcIO()); 2347 2348 using dof_map_t = std::map<types::global_dof_index, Point<spacedim>>; 2349 2350 using point_map_t = std::map<Point<spacedim>, 2351 std::vector<types::global_dof_index>, 2352 typename internal::ComparisonHelper<spacedim>>; 2353 2354 point_map_t point_map; 2355 2356 // convert to map point -> list of DoFs 2357 for (typename dof_map_t::const_iterator it = support_points.begin(); 2358 it != support_points.end(); 2359 ++it) 2360 { 2361 std::vector<types::global_dof_index> &v = point_map[it->second]; 2362 v.push_back(it->first); 2363 } 2364 2365 // print the newly created map: 2366 for (typename point_map_t::iterator it = point_map.begin(); 2367 it != point_map.end(); 2368 ++it) 2369 { 2370 out << it->first << " \""; 2371 const std::vector<types::global_dof_index> &v = it->second; 2372 for (unsigned int i = 0; i < v.size(); ++i) 2373 { 2374 if (i > 0) 2375 out << ", "; 2376 out << v[i]; 2377 } 2378 2379 out << "\"\n"; 2380 } 2381 2382 out << std::flush; 2383 } 2384 2385 2386 template <int dim, int spacedim> 2387 void convert_couplings_to_blocks(const DoFHandler<dim,spacedim> & dof_handler,const Table<2,Coupling> & table,std::vector<Table<2,Coupling>> & tables_by_block)2388 convert_couplings_to_blocks(const DoFHandler<dim, spacedim> &dof_handler, 2389 const Table<2, Coupling> & table, 2390 std::vector<Table<2, Coupling>> &tables_by_block) 2391 { 2392 if (dof_handler.hp_capability_enabled == false) 2393 { 2394 const FiniteElement<dim, spacedim> &fe = dof_handler.get_fe(); 2395 const unsigned int nb = fe.n_blocks(); 2396 2397 tables_by_block.resize(1); 2398 tables_by_block[0].reinit(nb, nb); 2399 tables_by_block[0].fill(none); 2400 2401 for (unsigned int i = 0; i < fe.n_components(); ++i) 2402 { 2403 const unsigned int ib = fe.component_to_block_index(i); 2404 for (unsigned int j = 0; j < fe.n_components(); ++j) 2405 { 2406 const unsigned int jb = fe.component_to_block_index(j); 2407 tables_by_block[0](ib, jb) |= table(i, j); 2408 } 2409 } 2410 } 2411 else 2412 { 2413 const hp::FECollection<dim> &fe_collection = 2414 dof_handler.get_fe_collection(); 2415 tables_by_block.resize(fe_collection.size()); 2416 2417 for (unsigned int f = 0; f < fe_collection.size(); ++f) 2418 { 2419 const FiniteElement<dim, spacedim> &fe = fe_collection[f]; 2420 2421 const unsigned int nb = fe.n_blocks(); 2422 tables_by_block[f].reinit(nb, nb); 2423 tables_by_block[f].fill(none); 2424 for (unsigned int i = 0; i < fe.n_components(); ++i) 2425 { 2426 const unsigned int ib = fe.component_to_block_index(i); 2427 for (unsigned int j = 0; j < fe.n_components(); ++j) 2428 { 2429 const unsigned int jb = fe.component_to_block_index(j); 2430 tables_by_block[f](ib, jb) |= table(i, j); 2431 } 2432 } 2433 } 2434 } 2435 } 2436 2437 2438 2439 template <int dim, int spacedim> 2440 void make_cell_patches(SparsityPattern & block_list,const DoFHandler<dim,spacedim> & dof_handler,const unsigned int level,const std::vector<bool> & selected_dofs,const types::global_dof_index offset)2441 make_cell_patches(SparsityPattern & block_list, 2442 const DoFHandler<dim, spacedim> &dof_handler, 2443 const unsigned int level, 2444 const std::vector<bool> & selected_dofs, 2445 const types::global_dof_index offset) 2446 { 2447 typename DoFHandler<dim, spacedim>::level_cell_iterator cell; 2448 typename DoFHandler<dim, spacedim>::level_cell_iterator endc = 2449 dof_handler.end(level); 2450 std::vector<types::global_dof_index> indices; 2451 2452 unsigned int i = 0; 2453 2454 for (cell = dof_handler.begin(level); cell != endc; ++cell) 2455 if (cell->is_locally_owned_on_level()) 2456 ++i; 2457 block_list.reinit(i, 2458 dof_handler.n_dofs(), 2459 dof_handler.get_fe().n_dofs_per_cell()); 2460 i = 0; 2461 for (cell = dof_handler.begin(level); cell != endc; ++cell) 2462 if (cell->is_locally_owned_on_level()) 2463 { 2464 indices.resize(cell->get_fe().n_dofs_per_cell()); 2465 cell->get_mg_dof_indices(indices); 2466 2467 if (selected_dofs.size() != 0) 2468 AssertDimension(indices.size(), selected_dofs.size()); 2469 2470 for (types::global_dof_index j = 0; j < indices.size(); ++j) 2471 { 2472 if (selected_dofs.size() == 0) 2473 block_list.add(i, indices[j] - offset); 2474 else 2475 { 2476 if (selected_dofs[j]) 2477 block_list.add(i, indices[j] - offset); 2478 } 2479 } 2480 ++i; 2481 } 2482 } 2483 2484 2485 template <int dim, int spacedim> 2486 void make_single_patch(SparsityPattern & block_list,const DoFHandler<dim,spacedim> & dof_handler,const unsigned int level,const bool interior_only)2487 make_single_patch(SparsityPattern & block_list, 2488 const DoFHandler<dim, spacedim> &dof_handler, 2489 const unsigned int level, 2490 const bool interior_only) 2491 { 2492 const FiniteElement<dim> &fe = dof_handler.get_fe(); 2493 block_list.reinit(1, dof_handler.n_dofs(level), dof_handler.n_dofs(level)); 2494 typename DoFHandler<dim, spacedim>::level_cell_iterator cell; 2495 typename DoFHandler<dim, spacedim>::level_cell_iterator endc = 2496 dof_handler.end(level); 2497 2498 std::vector<types::global_dof_index> indices; 2499 std::vector<bool> exclude; 2500 2501 for (cell = dof_handler.begin(level); cell != endc; ++cell) 2502 { 2503 indices.resize(cell->get_fe().n_dofs_per_cell()); 2504 cell->get_mg_dof_indices(indices); 2505 2506 if (interior_only) 2507 { 2508 // Exclude degrees of freedom on faces opposite to the vertex 2509 exclude.resize(fe.n_dofs_per_cell()); 2510 std::fill(exclude.begin(), exclude.end(), false); 2511 const unsigned int dpf = fe.n_dofs_per_face(); 2512 2513 for (const unsigned int face : cell->face_indices()) 2514 if (cell->at_boundary(face) || 2515 cell->neighbor(face)->level() != cell->level()) 2516 for (unsigned int i = 0; i < dpf; ++i) 2517 exclude[fe.face_to_cell_index(i, face)] = true; 2518 for (types::global_dof_index j = 0; j < indices.size(); ++j) 2519 if (!exclude[j]) 2520 block_list.add(0, indices[j]); 2521 } 2522 else 2523 { 2524 for (const auto index : indices) 2525 block_list.add(0, index); 2526 } 2527 } 2528 } 2529 2530 2531 template <int dim, int spacedim> 2532 void make_child_patches(SparsityPattern & block_list,const DoFHandler<dim,spacedim> & dof_handler,const unsigned int level,const bool interior_dofs_only,const bool boundary_dofs)2533 make_child_patches(SparsityPattern & block_list, 2534 const DoFHandler<dim, spacedim> &dof_handler, 2535 const unsigned int level, 2536 const bool interior_dofs_only, 2537 const bool boundary_dofs) 2538 { 2539 Assert(level > 0 && level < dof_handler.get_triangulation().n_levels(), 2540 ExcIndexRange(level, 1, dof_handler.get_triangulation().n_levels())); 2541 2542 typename DoFHandler<dim, spacedim>::level_cell_iterator pcell = 2543 dof_handler.begin(level - 1); 2544 typename DoFHandler<dim, spacedim>::level_cell_iterator endc = 2545 dof_handler.end(level - 1); 2546 2547 std::vector<types::global_dof_index> indices; 2548 std::vector<bool> exclude; 2549 2550 for (unsigned int block = 0; pcell != endc; ++pcell) 2551 { 2552 if (pcell->is_active()) 2553 continue; 2554 2555 for (unsigned int child = 0; child < pcell->n_children(); ++child) 2556 { 2557 const typename DoFHandler<dim, spacedim>::level_cell_iterator cell = 2558 pcell->child(child); 2559 2560 // For hp, only this line here would have to be replaced. 2561 const FiniteElement<dim> &fe = dof_handler.get_fe(); 2562 const unsigned int n_dofs = fe.n_dofs_per_cell(); 2563 indices.resize(n_dofs); 2564 exclude.resize(n_dofs); 2565 std::fill(exclude.begin(), exclude.end(), false); 2566 cell->get_mg_dof_indices(indices); 2567 2568 if (interior_dofs_only) 2569 { 2570 // Eliminate dofs on faces of the child which are on faces 2571 // of the parent 2572 const unsigned int dpf = fe.n_dofs_per_face(); 2573 2574 for (unsigned int d = 0; d < dim; ++d) 2575 { 2576 const unsigned int face = 2577 GeometryInfo<dim>::vertex_to_face[child][d]; 2578 for (unsigned int i = 0; i < dpf; ++i) 2579 exclude[fe.face_to_cell_index(i, face)] = true; 2580 } 2581 2582 // Now remove all degrees of freedom on the domain boundary 2583 // from the exclusion list 2584 if (boundary_dofs) 2585 for (const unsigned int face : 2586 GeometryInfo<dim>::face_indices()) 2587 if (cell->at_boundary(face)) 2588 for (unsigned int i = 0; i < dpf; ++i) 2589 exclude[fe.face_to_cell_index(i, face)] = false; 2590 } 2591 2592 for (unsigned int i = 0; i < n_dofs; ++i) 2593 if (!exclude[i]) 2594 block_list.add(block, indices[i]); 2595 } 2596 ++block; 2597 } 2598 } 2599 2600 template <int dim, int spacedim> 2601 std::vector<unsigned int> make_vertex_patches(SparsityPattern & block_list,const DoFHandler<dim,spacedim> & dof_handler,const unsigned int level,const bool interior_only,const bool boundary_patches,const bool level_boundary_patches,const bool single_cell_patches,const bool invert_vertex_mapping)2602 make_vertex_patches(SparsityPattern & block_list, 2603 const DoFHandler<dim, spacedim> &dof_handler, 2604 const unsigned int level, 2605 const bool interior_only, 2606 const bool boundary_patches, 2607 const bool level_boundary_patches, 2608 const bool single_cell_patches, 2609 const bool invert_vertex_mapping) 2610 { 2611 const unsigned int n_blocks = dof_handler.get_fe().n_blocks(); 2612 BlockMask exclude_boundary_dofs = BlockMask(n_blocks, interior_only); 2613 return make_vertex_patches(block_list, 2614 dof_handler, 2615 level, 2616 exclude_boundary_dofs, 2617 boundary_patches, 2618 level_boundary_patches, 2619 single_cell_patches, 2620 invert_vertex_mapping); 2621 } 2622 2623 template <int dim, int spacedim> 2624 std::vector<unsigned int> make_vertex_patches(SparsityPattern & block_list,const DoFHandler<dim,spacedim> & dof_handler,const unsigned int level,const BlockMask & exclude_boundary_dofs,const bool boundary_patches,const bool level_boundary_patches,const bool single_cell_patches,const bool invert_vertex_mapping)2625 make_vertex_patches(SparsityPattern & block_list, 2626 const DoFHandler<dim, spacedim> &dof_handler, 2627 const unsigned int level, 2628 const BlockMask & exclude_boundary_dofs, 2629 const bool boundary_patches, 2630 const bool level_boundary_patches, 2631 const bool single_cell_patches, 2632 const bool invert_vertex_mapping) 2633 { 2634 typename DoFHandler<dim, spacedim>::level_cell_iterator cell; 2635 typename DoFHandler<dim, spacedim>::level_cell_iterator endc = 2636 dof_handler.end(level); 2637 2638 // Vector mapping from vertex index in the triangulation to consecutive 2639 // block indices on this level The number of cells at a vertex 2640 std::vector<unsigned int> vertex_cell_count( 2641 dof_handler.get_triangulation().n_vertices(), 0); 2642 2643 // Is a vertex at the boundary? 2644 std::vector<bool> vertex_boundary( 2645 dof_handler.get_triangulation().n_vertices(), false); 2646 2647 std::vector<unsigned int> vertex_mapping( 2648 dof_handler.get_triangulation().n_vertices(), 2649 numbers::invalid_unsigned_int); 2650 2651 // Estimate for the number of dofs at this point 2652 std::vector<unsigned int> vertex_dof_count( 2653 dof_handler.get_triangulation().n_vertices(), 0); 2654 2655 // Identify all vertices active on this level and remember some data 2656 // about them 2657 for (cell = dof_handler.begin(level); cell != endc; ++cell) 2658 for (const unsigned int v : cell->vertex_indices()) 2659 { 2660 const unsigned int vg = cell->vertex_index(v); 2661 vertex_dof_count[vg] += cell->get_fe().n_dofs_per_cell(); 2662 ++vertex_cell_count[vg]; 2663 for (unsigned int d = 0; d < dim; ++d) 2664 { 2665 const unsigned int face = GeometryInfo<dim>::vertex_to_face[v][d]; 2666 if (cell->at_boundary(face)) 2667 vertex_boundary[vg] = true; 2668 else if ((!level_boundary_patches) && 2669 (cell->neighbor(face)->level() != 2670 static_cast<int>(level))) 2671 vertex_boundary[vg] = true; 2672 } 2673 } 2674 // From now on, only vertices with positive dof count are "in". 2675 2676 // Remove vertices at boundaries or in corners 2677 for (unsigned int vg = 0; vg < vertex_dof_count.size(); ++vg) 2678 if ((!single_cell_patches && vertex_cell_count[vg] < 2) || 2679 (!boundary_patches && vertex_boundary[vg])) 2680 vertex_dof_count[vg] = 0; 2681 2682 // Create a mapping from all vertices to the ones used here 2683 unsigned int n_vertex_count = 0; 2684 for (unsigned int vg = 0; vg < vertex_mapping.size(); ++vg) 2685 if (vertex_dof_count[vg] != 0) 2686 vertex_mapping[vg] = n_vertex_count++; 2687 2688 // Compactify dof count 2689 for (unsigned int vg = 0; vg < vertex_mapping.size(); ++vg) 2690 if (vertex_dof_count[vg] != 0) 2691 vertex_dof_count[vertex_mapping[vg]] = vertex_dof_count[vg]; 2692 2693 // Now that we have all the data, we reduce it to the part we actually 2694 // want 2695 vertex_dof_count.resize(n_vertex_count); 2696 2697 // At this point, the list of patches is ready. Now we enter the dofs 2698 // into the sparsity pattern. 2699 block_list.reinit(vertex_dof_count.size(), 2700 dof_handler.n_dofs(level), 2701 vertex_dof_count); 2702 2703 std::vector<types::global_dof_index> indices; 2704 std::vector<bool> exclude; 2705 2706 for (cell = dof_handler.begin(level); cell != endc; ++cell) 2707 { 2708 const FiniteElement<dim> &fe = cell->get_fe(); 2709 indices.resize(fe.n_dofs_per_cell()); 2710 cell->get_mg_dof_indices(indices); 2711 2712 for (const unsigned int v : cell->vertex_indices()) 2713 { 2714 const unsigned int vg = cell->vertex_index(v); 2715 const unsigned int block = vertex_mapping[vg]; 2716 if (block == numbers::invalid_unsigned_int) 2717 continue; 2718 2719 // Collect excluded dofs for some block(s) if boundary dofs 2720 // for a block are decided to be excluded 2721 if (exclude_boundary_dofs.size() == 0 || 2722 exclude_boundary_dofs.n_selected_blocks() != 0) 2723 { 2724 // Exclude degrees of freedom on faces opposite to the 2725 // vertex 2726 exclude.resize(fe.n_dofs_per_cell()); 2727 std::fill(exclude.begin(), exclude.end(), false); 2728 const unsigned int dpf = fe.n_dofs_per_face(); 2729 2730 for (unsigned int d = 0; d < dim; ++d) 2731 { 2732 const unsigned int a_face = 2733 GeometryInfo<dim>::vertex_to_face[v][d]; 2734 const unsigned int face = 2735 GeometryInfo<dim>::opposite_face[a_face]; 2736 for (unsigned int i = 0; i < dpf; ++i) 2737 { 2738 // For each dof, get the block it is in and decide to 2739 // exclude it or not 2740 if (exclude_boundary_dofs[fe.system_to_block_index( 2741 fe.face_to_cell_index( 2742 i, face)) 2743 .first] == true) 2744 exclude[fe.face_to_cell_index(i, face)] = true; 2745 } 2746 } 2747 for (unsigned int j = 0; j < indices.size(); ++j) 2748 if (!exclude[j]) 2749 block_list.add(block, indices[j]); 2750 } 2751 else 2752 { 2753 for (const auto index : indices) 2754 block_list.add(block, index); 2755 } 2756 } 2757 } 2758 2759 if (invert_vertex_mapping) 2760 { 2761 // Compress vertex mapping 2762 unsigned int n_vertex_count = 0; 2763 for (unsigned int vg = 0; vg < vertex_mapping.size(); ++vg) 2764 if (vertex_mapping[vg] != numbers::invalid_unsigned_int) 2765 vertex_mapping[n_vertex_count++] = vg; 2766 2767 // Now we reduce it to the part we actually want 2768 vertex_mapping.resize(n_vertex_count); 2769 } 2770 2771 return vertex_mapping; 2772 } 2773 2774 2775 template <int dim, int spacedim> 2776 unsigned int count_dofs_on_patch(const std::vector<typename DoFHandler<dim,spacedim>::active_cell_iterator> & patch)2777 count_dofs_on_patch( 2778 const std::vector<typename DoFHandler<dim, spacedim>::active_cell_iterator> 2779 &patch) 2780 { 2781 std::set<types::global_dof_index> dofs_on_patch; 2782 std::vector<types::global_dof_index> local_dof_indices; 2783 2784 // loop over the cells in the patch and get the DoFs on each. 2785 // add all of them to a std::set which automatically makes sure 2786 // all duplicates are ignored 2787 for (unsigned int i = 0; i < patch.size(); ++i) 2788 { 2789 const typename DoFHandler<dim, spacedim>::active_cell_iterator cell = 2790 patch[i]; 2791 Assert(cell->is_artificial() == false, 2792 ExcMessage("This function can not be called with cells that are " 2793 "not either locally owned or ghost cells.")); 2794 local_dof_indices.resize(cell->get_fe().n_dofs_per_cell()); 2795 cell->get_dof_indices(local_dof_indices); 2796 dofs_on_patch.insert(local_dof_indices.begin(), 2797 local_dof_indices.end()); 2798 } 2799 2800 // now return the number of DoFs (duplicates were ignored) 2801 return dofs_on_patch.size(); 2802 } 2803 2804 2805 2806 template <int dim, int spacedim> 2807 std::vector<types::global_dof_index> get_dofs_on_patch(const std::vector<typename DoFHandler<dim,spacedim>::active_cell_iterator> & patch)2808 get_dofs_on_patch( 2809 const std::vector<typename DoFHandler<dim, spacedim>::active_cell_iterator> 2810 &patch) 2811 { 2812 std::set<types::global_dof_index> dofs_on_patch; 2813 std::vector<types::global_dof_index> local_dof_indices; 2814 2815 // loop over the cells in the patch and get the DoFs on each. 2816 // add all of them to a std::set which automatically makes sure 2817 // all duplicates are ignored 2818 for (unsigned int i = 0; i < patch.size(); ++i) 2819 { 2820 const typename DoFHandler<dim, spacedim>::active_cell_iterator cell = 2821 patch[i]; 2822 Assert(cell->is_artificial() == false, 2823 ExcMessage("This function can not be called with cells that are " 2824 "not either locally owned or ghost cells.")); 2825 local_dof_indices.resize(cell->get_fe().n_dofs_per_cell()); 2826 cell->get_dof_indices(local_dof_indices); 2827 dofs_on_patch.insert(local_dof_indices.begin(), 2828 local_dof_indices.end()); 2829 } 2830 2831 Assert((dofs_on_patch.size() == count_dofs_on_patch<dim, spacedim>(patch)), 2832 ExcInternalError()); 2833 2834 // return a vector with the content of the set above. copying 2835 // also ensures that we retain sortedness as promised in the 2836 // documentation and as necessary to retain the block structure 2837 // also on the local system 2838 return std::vector<types::global_dof_index>(dofs_on_patch.begin(), 2839 dofs_on_patch.end()); 2840 } 2841 2842 2843 } // end of namespace DoFTools 2844 2845 2846 2847 // explicit instantiations 2848 2849 #include "dof_tools.inst" 2850 2851 2852 2853 DEAL_II_NAMESPACE_CLOSE 2854