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/table.h> 17 #include <deal.II/base/template_constraints.h> 18 #include <deal.II/base/utilities.h> 19 #include <deal.II/base/work_stream.h> 20 21 #include <deal.II/dofs/dof_accessor.h> 22 #include <deal.II/dofs/dof_handler.h> 23 #include <deal.II/dofs/dof_tools.h> 24 25 #include <deal.II/fe/fe.h> 26 #include <deal.II/fe/fe_tools.h> 27 #include <deal.II/fe/fe_values.h> 28 29 #include <deal.II/grid/grid_tools.h> 30 #include <deal.II/grid/intergrid_map.h> 31 #include <deal.II/grid/tria.h> 32 #include <deal.II/grid/tria_iterator.h> 33 34 #include <deal.II/hp/fe_collection.h> 35 #include <deal.II/hp/fe_values.h> 36 37 #include <deal.II/lac/affine_constraints.h> 38 #include <deal.II/lac/vector.h> 39 40 #ifdef DEAL_II_WITH_MPI 41 # include <deal.II/lac/la_parallel_vector.h> 42 #endif 43 44 #include <algorithm> 45 #include <array> 46 #include <memory> 47 #include <numeric> 48 49 DEAL_II_NAMESPACE_OPEN 50 51 52 53 namespace DoFTools 54 { 55 namespace internal 56 { 57 namespace 58 { 59 inline bool check_primary_dof_list(const FullMatrix<double> & face_interpolation_matrix,const std::vector<types::global_dof_index> & primary_dof_list)60 check_primary_dof_list( 61 const FullMatrix<double> & face_interpolation_matrix, 62 const std::vector<types::global_dof_index> &primary_dof_list) 63 { 64 const unsigned int N = primary_dof_list.size(); 65 66 FullMatrix<double> tmp(N, N); 67 for (unsigned int i = 0; i < N; ++i) 68 for (unsigned int j = 0; j < N; ++j) 69 tmp(i, j) = face_interpolation_matrix(primary_dof_list[i], j); 70 71 // then use the algorithm from FullMatrix::gauss_jordan on this matrix 72 // to find out whether it is singular. the algorithm there does pivoting 73 // and at the end swaps rows back into their proper order -- we omit 74 // this step here, since we don't care about the inverse matrix, all we 75 // care about is whether the matrix is regular or singular 76 77 // first get an estimate of the size of the elements of this matrix, for 78 // later checks whether the pivot element is large enough, or whether we 79 // have to fear that the matrix is not regular 80 double diagonal_sum = 0; 81 for (unsigned int i = 0; i < N; ++i) 82 diagonal_sum += std::fabs(tmp(i, i)); 83 const double typical_diagonal_element = diagonal_sum / N; 84 85 // initialize the array that holds the permutations that we find during 86 // pivot search 87 std::vector<unsigned int> p(N); 88 for (unsigned int i = 0; i < N; ++i) 89 p[i] = i; 90 91 for (unsigned int j = 0; j < N; ++j) 92 { 93 // pivot search: search that part of the line on and right of the 94 // diagonal for the largest element 95 double max = std::fabs(tmp(j, j)); 96 unsigned int r = j; 97 for (unsigned int i = j + 1; i < N; ++i) 98 { 99 if (std::fabs(tmp(i, j)) > max) 100 { 101 max = std::fabs(tmp(i, j)); 102 r = i; 103 } 104 } 105 // check whether the pivot is too small. if that is the case, then 106 // the matrix is singular and we shouldn't use this set of primary 107 // dofs 108 if (max < 1.e-12 * typical_diagonal_element) 109 return false; 110 111 // row interchange 112 if (r > j) 113 { 114 for (unsigned int k = 0; k < N; ++k) 115 std::swap(tmp(j, k), tmp(r, k)); 116 117 std::swap(p[j], p[r]); 118 } 119 120 // transformation 121 const double hr = 1. / tmp(j, j); 122 tmp(j, j) = hr; 123 for (unsigned int k = 0; k < N; ++k) 124 { 125 if (k == j) 126 continue; 127 for (unsigned int i = 0; i < N; ++i) 128 { 129 if (i == j) 130 continue; 131 tmp(i, k) -= tmp(i, j) * tmp(j, k) * hr; 132 } 133 } 134 for (unsigned int i = 0; i < N; ++i) 135 { 136 tmp(i, j) *= hr; 137 tmp(j, i) *= -hr; 138 } 139 tmp(j, j) = hr; 140 } 141 142 // everything went fine, so we can accept this set of primary dofs (at 143 // least as far as they have already been collected) 144 return true; 145 } 146 147 148 149 /** 150 * When restricting, on a face, the degrees of freedom of fe1 to the space 151 * described by fe2 (for example for the complex case described 152 * in the @ref hp_paper "hp paper"), we have to select fe2.dofs_per_face 153 * out of the fe1.dofs_per_face face DoFs as the 154 * primary dofs, and the rest become dependent dofs. This function selects 155 * which ones will be primary, and which ones will be dependents. 156 * 157 * The function assumes that primary_dofs already has size 158 * fe1.dofs_per_face. After the function, exactly fe2.dofs_per_face 159 * entries will be true. 160 * 161 * The function is a bit complicated since it has to figure out a set a 162 * DoFs so that the corresponding rows in the face interpolation matrix 163 * are all linearly independent. we have a good heuristic (see the 164 * function body) for selecting these rows, but there are cases where this 165 * fails and we have to pick them differently. what we do is to run the 166 * heuristic and then go back to determine whether we have a set of rows 167 * with full row rank. if this isn't the case, go back and select dofs 168 * differently 169 */ 170 template <int dim, int spacedim> 171 void select_primary_dofs_for_face_restriction(const FiniteElement<dim,spacedim> & fe1,const FiniteElement<dim,spacedim> & fe2,const FullMatrix<double> & face_interpolation_matrix,std::vector<bool> & primary_dof_mask)172 select_primary_dofs_for_face_restriction( 173 const FiniteElement<dim, spacedim> &fe1, 174 const FiniteElement<dim, spacedim> &fe2, 175 const FullMatrix<double> & face_interpolation_matrix, 176 std::vector<bool> & primary_dof_mask) 177 { 178 Assert(fe1.n_dofs_per_face() >= fe2.n_dofs_per_face(), 179 ExcInternalError()); 180 AssertDimension(primary_dof_mask.size(), fe1.n_dofs_per_face()); 181 182 Assert(fe2.n_dofs_per_vertex() <= fe1.n_dofs_per_vertex(), 183 ExcInternalError()); 184 Assert(fe2.n_dofs_per_line() <= fe1.n_dofs_per_line(), 185 ExcInternalError()); 186 Assert((dim < 3) || (fe2.n_dofs_per_quad() <= fe1.n_dofs_per_quad()), 187 ExcInternalError()); 188 189 // the idea here is to designate as many DoFs in fe1 per object (vertex, 190 // line, quad) as primary as there are such dofs in fe2 (indices are 191 // int, because we want to avoid the 'unsigned int < 0 is always false 192 // warning for the cases at the bottom in 1d and 2d) 193 // 194 // as mentioned in the paper, it is not always easy to find a set of 195 // primary dofs that produces an invertible matrix. to this end, we 196 // check in each step whether the matrix is still invertible and simply 197 // discard this dof if the matrix is not invertible anymore. 198 // 199 // the cases where we did have trouble in the past were with adding more 200 // quad dofs when Q3 and Q4 elements meet at a refined face in 3d (see 201 // the hp/crash_12 test that tests that we can do exactly this, and 202 // failed before we had code to compensate for this case). the other 203 // case are system elements: if we have say a Q1Q2 vs a Q2Q3 element, 204 // then we can't just take all primary dofs on a line from a single base 205 // element, since the shape functions of that base element are 206 // independent of that of the other one. this latter case shows up when 207 // running hp/hp_constraints_q_system_06 208 209 std::vector<types::global_dof_index> primary_dof_list; 210 unsigned int index = 0; 211 for (int v = 0; 212 v < static_cast<signed int>(GeometryInfo<dim>::vertices_per_face); 213 ++v) 214 { 215 unsigned int dofs_added = 0; 216 unsigned int i = 0; 217 while (dofs_added < fe2.n_dofs_per_vertex()) 218 { 219 // make sure that we were able to find a set of primary dofs and 220 // that the code down below didn't just reject all our efforts 221 Assert(i < fe1.n_dofs_per_vertex(), ExcInternalError()); 222 223 // tentatively push this vertex dof 224 primary_dof_list.push_back(index + i); 225 226 // then see what happens. if it succeeds, fine 227 if (check_primary_dof_list(face_interpolation_matrix, 228 primary_dof_list) == true) 229 ++dofs_added; 230 else 231 // well, it didn't. simply pop that dof from the list again 232 // and try with the next dof 233 primary_dof_list.pop_back(); 234 235 // forward counter by one 236 ++i; 237 } 238 index += fe1.n_dofs_per_vertex(); 239 } 240 241 for (int l = 0; 242 l < static_cast<signed int>(GeometryInfo<dim>::lines_per_face); 243 ++l) 244 { 245 // same algorithm as above 246 unsigned int dofs_added = 0; 247 unsigned int i = 0; 248 while (dofs_added < fe2.n_dofs_per_line()) 249 { 250 Assert(i < fe1.n_dofs_per_line(), ExcInternalError()); 251 252 primary_dof_list.push_back(index + i); 253 if (check_primary_dof_list(face_interpolation_matrix, 254 primary_dof_list) == true) 255 ++dofs_added; 256 else 257 primary_dof_list.pop_back(); 258 259 ++i; 260 } 261 index += fe1.n_dofs_per_line(); 262 } 263 264 for (int q = 0; 265 q < static_cast<signed int>(GeometryInfo<dim>::quads_per_face); 266 ++q) 267 { 268 // same algorithm as above 269 unsigned int dofs_added = 0; 270 unsigned int i = 0; 271 while (dofs_added < fe2.n_dofs_per_quad()) 272 { 273 Assert(i < fe1.n_dofs_per_quad(), ExcInternalError()); 274 275 primary_dof_list.push_back(index + i); 276 if (check_primary_dof_list(face_interpolation_matrix, 277 primary_dof_list) == true) 278 ++dofs_added; 279 else 280 primary_dof_list.pop_back(); 281 282 ++i; 283 } 284 index += fe1.n_dofs_per_quad(); 285 } 286 287 AssertDimension(index, fe1.n_dofs_per_face()); 288 AssertDimension(primary_dof_list.size(), fe2.n_dofs_per_face()); 289 290 // finally copy the list into the mask 291 std::fill(primary_dof_mask.begin(), primary_dof_mask.end(), false); 292 for (const auto dof : primary_dof_list) 293 primary_dof_mask[dof] = true; 294 } 295 296 297 298 /** 299 * Make sure that the mask exists that determines which dofs will be the 300 * primary on refined faces where an fe1 and a fe2 meet. 301 */ 302 template <int dim, int spacedim> 303 void ensure_existence_of_primary_dof_mask(const FiniteElement<dim,spacedim> & fe1,const FiniteElement<dim,spacedim> & fe2,const FullMatrix<double> & face_interpolation_matrix,std::unique_ptr<std::vector<bool>> & primary_dof_mask)304 ensure_existence_of_primary_dof_mask( 305 const FiniteElement<dim, spacedim> &fe1, 306 const FiniteElement<dim, spacedim> &fe2, 307 const FullMatrix<double> & face_interpolation_matrix, 308 std::unique_ptr<std::vector<bool>> &primary_dof_mask) 309 { 310 if (primary_dof_mask == nullptr) 311 { 312 primary_dof_mask = 313 std::make_unique<std::vector<bool>>(fe1.n_dofs_per_face()); 314 select_primary_dofs_for_face_restriction(fe1, 315 fe2, 316 face_interpolation_matrix, 317 *primary_dof_mask); 318 } 319 } 320 321 322 323 /** 324 * Make sure that the given @p face_interpolation_matrix pointer points 325 * to a valid matrix. If the pointer is zero beforehand, create an entry 326 * with the correct data. If it is nonzero, don't touch it. 327 */ 328 template <int dim, int spacedim> 329 void ensure_existence_of_face_matrix(const FiniteElement<dim,spacedim> & fe1,const FiniteElement<dim,spacedim> & fe2,std::unique_ptr<FullMatrix<double>> & matrix)330 ensure_existence_of_face_matrix( 331 const FiniteElement<dim, spacedim> & fe1, 332 const FiniteElement<dim, spacedim> & fe2, 333 std::unique_ptr<FullMatrix<double>> &matrix) 334 { 335 if (matrix == nullptr) 336 { 337 matrix = 338 std::make_unique<FullMatrix<double>>(fe2.n_dofs_per_face(), 339 fe1.n_dofs_per_face()); 340 fe1.get_face_interpolation_matrix(fe2, *matrix); 341 } 342 } 343 344 345 346 /** 347 * Same, but for subface interpolation matrices. 348 */ 349 template <int dim, int spacedim> 350 void ensure_existence_of_subface_matrix(const FiniteElement<dim,spacedim> & fe1,const FiniteElement<dim,spacedim> & fe2,const unsigned int subface,std::unique_ptr<FullMatrix<double>> & matrix)351 ensure_existence_of_subface_matrix( 352 const FiniteElement<dim, spacedim> & fe1, 353 const FiniteElement<dim, spacedim> & fe2, 354 const unsigned int subface, 355 std::unique_ptr<FullMatrix<double>> &matrix) 356 { 357 if (matrix == nullptr) 358 { 359 matrix = 360 std::make_unique<FullMatrix<double>>(fe2.n_dofs_per_face(), 361 fe1.n_dofs_per_face()); 362 fe1.get_subface_interpolation_matrix(fe2, subface, *matrix); 363 } 364 } 365 366 367 368 /** 369 * Given the face interpolation matrix between two elements, split it into 370 * its primary and dependent parts and invert the primary part as 371 * explained in the @ref hp_paper "hp paper". 372 */ 373 void ensure_existence_of_split_face_matrix(const FullMatrix<double> & face_interpolation_matrix,const std::vector<bool> & primary_dof_mask,std::unique_ptr<std::pair<FullMatrix<double>,FullMatrix<double>>> & split_matrix)374 ensure_existence_of_split_face_matrix( 375 const FullMatrix<double> &face_interpolation_matrix, 376 const std::vector<bool> & primary_dof_mask, 377 std::unique_ptr<std::pair<FullMatrix<double>, FullMatrix<double>>> 378 &split_matrix) 379 { 380 AssertDimension(primary_dof_mask.size(), face_interpolation_matrix.m()); 381 Assert(std::count(primary_dof_mask.begin(), 382 primary_dof_mask.end(), 383 true) == 384 static_cast<signed int>(face_interpolation_matrix.n()), 385 ExcInternalError()); 386 387 if (split_matrix == nullptr) 388 { 389 split_matrix = std::make_unique< 390 std::pair<FullMatrix<double>, FullMatrix<double>>>(); 391 392 const unsigned int n_primary_dofs = face_interpolation_matrix.n(); 393 const unsigned int n_dofs = face_interpolation_matrix.m(); 394 395 Assert(n_primary_dofs <= n_dofs, ExcInternalError()); 396 397 // copy and invert the primary component, copy the dependent 398 // component 399 split_matrix->first.reinit(n_primary_dofs, n_primary_dofs); 400 split_matrix->second.reinit(n_dofs - n_primary_dofs, 401 n_primary_dofs); 402 403 unsigned int nth_primary_dof = 0, nth_dependent_dof = 0; 404 405 for (unsigned int i = 0; i < n_dofs; ++i) 406 if (primary_dof_mask[i] == true) 407 { 408 for (unsigned int j = 0; j < n_primary_dofs; ++j) 409 split_matrix->first(nth_primary_dof, j) = 410 face_interpolation_matrix(i, j); 411 ++nth_primary_dof; 412 } 413 else 414 { 415 for (unsigned int j = 0; j < n_primary_dofs; ++j) 416 split_matrix->second(nth_dependent_dof, j) = 417 face_interpolation_matrix(i, j); 418 ++nth_dependent_dof; 419 } 420 421 AssertDimension(nth_primary_dof, n_primary_dofs); 422 AssertDimension(nth_dependent_dof, n_dofs - n_primary_dofs); 423 424 // TODO[WB]: We should make sure very small entries are removed 425 // after inversion 426 split_matrix->first.gauss_jordan(); 427 } 428 } 429 430 431 /** 432 * A function that returns how many different finite elements a dof 433 * handler uses. This is one for non-hp DoFHandlers and 434 * dof_handler.get_fe().size() for the hp-versions. 435 */ 436 template <int dim, int spacedim> 437 unsigned int n_finite_elements(const DoFHandler<dim,spacedim> & dof_handler)438 n_finite_elements(const DoFHandler<dim, spacedim> &dof_handler) 439 { 440 if (dof_handler.hp_capability_enabled == true) 441 return dof_handler.get_fe_collection().size(); 442 else 443 return 1; 444 } 445 446 447 448 /** 449 * Copy constraints into an AffineConstraints object. 450 * 451 * This function removes zero constraints and those, which constrain a DoF 452 * which was already eliminated in one of the previous steps of the hp 453 * hanging node procedure. 454 * 455 * It also suppresses very small entries in the AffineConstraints object 456 * to avoid making the sparsity pattern fuller than necessary. 457 */ 458 template <typename number1, typename number2> 459 void filter_constraints(const std::vector<types::global_dof_index> & primary_dofs,const std::vector<types::global_dof_index> & dependent_dofs,const FullMatrix<number1> & face_constraints,AffineConstraints<number2> & constraints)460 filter_constraints( 461 const std::vector<types::global_dof_index> &primary_dofs, 462 const std::vector<types::global_dof_index> &dependent_dofs, 463 const FullMatrix<number1> & face_constraints, 464 AffineConstraints<number2> & constraints) 465 { 466 Assert(face_constraints.n() == primary_dofs.size(), 467 ExcDimensionMismatch(primary_dofs.size(), face_constraints.n())); 468 Assert(face_constraints.m() == dependent_dofs.size(), 469 ExcDimensionMismatch(dependent_dofs.size(), 470 face_constraints.m())); 471 472 const unsigned int n_primary_dofs = primary_dofs.size(); 473 const unsigned int n_dependent_dofs = dependent_dofs.size(); 474 475 // check for a couple conditions that happened in parallel distributed 476 // mode 477 for (unsigned int row = 0; row != n_dependent_dofs; ++row) 478 Assert(dependent_dofs[row] != numbers::invalid_dof_index, 479 ExcInternalError()); 480 for (unsigned int col = 0; col != n_primary_dofs; ++col) 481 Assert(primary_dofs[col] != numbers::invalid_dof_index, 482 ExcInternalError()); 483 484 485 for (unsigned int row = 0; row != n_dependent_dofs; ++row) 486 if (constraints.is_constrained(dependent_dofs[row]) == false) 487 { 488 bool constraint_already_satisfied = false; 489 490 // Check if we have an identity constraint, which is already 491 // satisfied by unification of the corresponding global dof 492 // indices 493 for (unsigned int i = 0; i < n_primary_dofs; ++i) 494 if (face_constraints(row, i) == 1.0) 495 if (primary_dofs[i] == dependent_dofs[row]) 496 { 497 constraint_already_satisfied = true; 498 break; 499 } 500 501 if (constraint_already_satisfied == false) 502 { 503 // add up the absolute values of all constraints in this line 504 // to get a measure of their absolute size 505 number1 abs_sum = 0; 506 for (unsigned int i = 0; i < n_primary_dofs; ++i) 507 abs_sum += std::abs(face_constraints(row, i)); 508 509 // then enter those constraints that are larger than 510 // 1e-14*abs_sum. everything else probably originated from 511 // inexact inversion of matrices and similar effects. having 512 // those constraints in here will only lead to problems 513 // because it makes sparsity patterns fuller than necessary 514 // without producing any significant effect 515 constraints.add_line(dependent_dofs[row]); 516 for (unsigned int i = 0; i < n_primary_dofs; ++i) 517 if ((face_constraints(row, i) != 0) && 518 (std::fabs(face_constraints(row, i)) >= 519 1e-14 * abs_sum)) 520 constraints.add_entry(dependent_dofs[row], 521 primary_dofs[i], 522 face_constraints(row, i)); 523 constraints.set_inhomogeneity(dependent_dofs[row], 0.); 524 } 525 } 526 } 527 528 } // namespace 529 530 531 template <typename number> 532 void make_hp_hanging_node_constraints(const dealii::DoFHandler<1> &,AffineConstraints<number> &)533 make_hp_hanging_node_constraints(const dealii::DoFHandler<1> &, 534 AffineConstraints<number> &) 535 { 536 // nothing to do for regular dof handlers in 1d 537 } 538 539 540 template <typename number> 541 void make_oldstyle_hanging_node_constraints(const dealii::DoFHandler<1> &,AffineConstraints<number> &,std::integral_constant<int,1>)542 make_oldstyle_hanging_node_constraints(const dealii::DoFHandler<1> &, 543 AffineConstraints<number> &, 544 std::integral_constant<int, 1>) 545 { 546 // nothing to do for regular dof handlers in 1d 547 } 548 549 550 template <typename number> 551 void make_hp_hanging_node_constraints(const dealii::DoFHandler<1,2> &,AffineConstraints<number> &)552 make_hp_hanging_node_constraints(const dealii::DoFHandler<1, 2> &, 553 AffineConstraints<number> &) 554 { 555 // nothing to do for regular dof handlers in 1d 556 } 557 558 559 template <typename number> 560 void make_oldstyle_hanging_node_constraints(const dealii::DoFHandler<1,2> &,AffineConstraints<number> &,std::integral_constant<int,1>)561 make_oldstyle_hanging_node_constraints(const dealii::DoFHandler<1, 2> &, 562 AffineConstraints<number> &, 563 std::integral_constant<int, 1>) 564 { 565 // nothing to do for regular dof handlers in 1d 566 } 567 568 569 template <typename number, int spacedim> 570 void make_hp_hanging_node_constraints(const dealii::DoFHandler<1,spacedim> &,AffineConstraints<number> &)571 make_hp_hanging_node_constraints( 572 const dealii::DoFHandler<1, spacedim> & /*dof_handler*/, 573 AffineConstraints<number> & /*constraints*/) 574 { 575 // nothing to do for dof handlers in 1d 576 } 577 578 579 template <typename number, int spacedim> 580 void make_oldstyle_hanging_node_constraints(const dealii::DoFHandler<1,spacedim> &,AffineConstraints<number> &,std::integral_constant<int,1>)581 make_oldstyle_hanging_node_constraints( 582 const dealii::DoFHandler<1, spacedim> & /*dof_handler*/, 583 AffineConstraints<number> & /*constraints*/, 584 std::integral_constant<int, 1>) 585 { 586 // nothing to do for dof handlers in 1d 587 } 588 589 template <int dim_, int spacedim, typename number> 590 void make_oldstyle_hanging_node_constraints(const DoFHandler<dim_,spacedim> & dof_handler,AffineConstraints<number> & constraints,std::integral_constant<int,2>)591 make_oldstyle_hanging_node_constraints( 592 const DoFHandler<dim_, spacedim> &dof_handler, 593 AffineConstraints<number> & constraints, 594 std::integral_constant<int, 2>) 595 { 596 const unsigned int dim = 2; 597 598 std::vector<types::global_dof_index> dofs_on_mother; 599 std::vector<types::global_dof_index> dofs_on_children; 600 601 // loop over all lines; only on lines there can be constraints. We do so 602 // by looping over all active cells and checking whether any of the faces 603 // are refined which can only be from the neighboring cell because this 604 // one is active. In that case, the face is subject to constraints 605 // 606 // note that even though we may visit a face twice if the neighboring 607 // cells are equally refined, we can only visit each face with hanging 608 // nodes once 609 typename DoFHandler<dim_, spacedim>::active_cell_iterator 610 cell = dof_handler.begin_active(), 611 endc = dof_handler.end(); 612 for (; cell != endc; ++cell) 613 { 614 // artificial cells can at best neighbor ghost cells, but we're not 615 // interested in these interfaces 616 if (cell->is_artificial()) 617 continue; 618 619 for (const unsigned int face : cell->face_indices()) 620 if (cell->face(face)->has_children()) 621 { 622 // in any case, faces can have at most two active fe indices, 623 // but here the face can have only one (namely the same as that 624 // from the cell we're sitting on), and each of the children can 625 // have only one as well. check this 626 Assert(cell->face(face)->n_active_fe_indices() == 1, 627 ExcInternalError()); 628 Assert(cell->face(face)->fe_index_is_active( 629 cell->active_fe_index()) == true, 630 ExcInternalError()); 631 for (unsigned int c = 0; c < cell->face(face)->n_children(); 632 ++c) 633 if (!cell->neighbor_child_on_subface(face, c) 634 ->is_artificial()) 635 Assert(cell->face(face)->child(c)->n_active_fe_indices() == 636 1, 637 ExcInternalError()); 638 639 // right now, all that is implemented is the case that both 640 // sides use the same fe 641 for (unsigned int c = 0; c < cell->face(face)->n_children(); 642 ++c) 643 if (!cell->neighbor_child_on_subface(face, c) 644 ->is_artificial()) 645 Assert(cell->face(face)->child(c)->fe_index_is_active( 646 cell->active_fe_index()) == true, 647 ExcNotImplemented()); 648 649 // ok, start up the work 650 const FiniteElement<dim, spacedim> &fe = cell->get_fe(); 651 const unsigned int fe_index = cell->active_fe_index(); 652 653 const unsigned int n_dofs_on_mother = 654 2 * fe.n_dofs_per_vertex() + 655 fe.n_dofs_per_line(), 656 n_dofs_on_children = 657 fe.n_dofs_per_vertex() + 658 2 * fe.n_dofs_per_line(); 659 660 dofs_on_mother.resize(n_dofs_on_mother); 661 // we might not use all of those in case of artificial cells, so 662 // do not resize(), but reserve() and use push_back later. 663 dofs_on_children.clear(); 664 dofs_on_children.reserve(n_dofs_on_children); 665 666 Assert(n_dofs_on_mother == fe.constraints().n(), 667 ExcDimensionMismatch(n_dofs_on_mother, 668 fe.constraints().n())); 669 Assert(n_dofs_on_children == fe.constraints().m(), 670 ExcDimensionMismatch(n_dofs_on_children, 671 fe.constraints().m())); 672 673 const typename DoFHandler<dim_, spacedim>::line_iterator 674 this_face = cell->face(face); 675 676 // fill the dofs indices. Use same enumeration scheme as in 677 // @p{FiniteElement::constraints()} 678 unsigned int next_index = 0; 679 for (unsigned int vertex = 0; vertex < 2; ++vertex) 680 for (unsigned int dof = 0; dof != fe.n_dofs_per_vertex(); 681 ++dof) 682 dofs_on_mother[next_index++] = 683 this_face->vertex_dof_index(vertex, dof, fe_index); 684 for (unsigned int dof = 0; dof != fe.n_dofs_per_line(); ++dof) 685 dofs_on_mother[next_index++] = 686 this_face->dof_index(dof, fe_index); 687 AssertDimension(next_index, dofs_on_mother.size()); 688 689 for (unsigned int dof = 0; dof != fe.n_dofs_per_vertex(); ++dof) 690 dofs_on_children.push_back( 691 this_face->child(0)->vertex_dof_index(1, dof, fe_index)); 692 for (unsigned int child = 0; child < 2; ++child) 693 { 694 // skip artificial cells 695 if (cell->neighbor_child_on_subface(face, child) 696 ->is_artificial()) 697 continue; 698 for (unsigned int dof = 0; dof != fe.n_dofs_per_line(); 699 ++dof) 700 dofs_on_children.push_back( 701 this_face->child(child)->dof_index(dof, fe_index)); 702 } 703 // note: can get fewer DoFs when we have artificial cells 704 Assert(dofs_on_children.size() <= n_dofs_on_children, 705 ExcInternalError()); 706 707 // for each row in the AffineConstraints object for this line: 708 for (unsigned int row = 0; row != dofs_on_children.size(); 709 ++row) 710 { 711 constraints.add_line(dofs_on_children[row]); 712 for (unsigned int i = 0; i != dofs_on_mother.size(); ++i) 713 constraints.add_entry(dofs_on_children[row], 714 dofs_on_mother[i], 715 fe.constraints()(row, i)); 716 717 constraints.set_inhomogeneity(dofs_on_children[row], 0.); 718 } 719 } 720 else 721 { 722 // this face has no children, but it could still be that it is 723 // shared by two cells that use a different fe index. check a 724 // couple of things, but ignore the case that the neighbor is an 725 // artificial cell 726 if (!cell->at_boundary(face) && 727 !cell->neighbor(face)->is_artificial()) 728 { 729 Assert(cell->face(face)->n_active_fe_indices() == 1, 730 ExcNotImplemented()); 731 Assert(cell->face(face)->fe_index_is_active( 732 cell->active_fe_index()) == true, 733 ExcInternalError()); 734 } 735 } 736 } 737 } 738 739 740 741 template <int dim_, int spacedim, typename number> 742 void make_oldstyle_hanging_node_constraints(const DoFHandler<dim_,spacedim> & dof_handler,AffineConstraints<number> & constraints,std::integral_constant<int,3>)743 make_oldstyle_hanging_node_constraints( 744 const DoFHandler<dim_, spacedim> &dof_handler, 745 AffineConstraints<number> & constraints, 746 std::integral_constant<int, 3>) 747 { 748 const unsigned int dim = 3; 749 750 std::vector<types::global_dof_index> dofs_on_mother; 751 std::vector<types::global_dof_index> dofs_on_children; 752 753 // loop over all quads; only on quads there can be constraints. We do so 754 // by looping over all active cells and checking whether any of the faces 755 // are refined which can only be from the neighboring cell because this 756 // one is active. In that case, the face is subject to constraints 757 // 758 // note that even though we may visit a face twice if the neighboring 759 // cells are equally refined, we can only visit each face with hanging 760 // nodes once 761 typename DoFHandler<dim_, spacedim>::active_cell_iterator 762 cell = dof_handler.begin_active(), 763 endc = dof_handler.end(); 764 for (; cell != endc; ++cell) 765 { 766 // artificial cells can at best neighbor ghost cells, but we're not 767 // interested in these interfaces 768 if (cell->is_artificial()) 769 continue; 770 771 for (const unsigned int face : cell->face_indices()) 772 if (cell->face(face)->has_children()) 773 { 774 // first of all, make sure that we treat a case which is 775 // possible, i.e. either no dofs on the face at all or no 776 // anisotropic refinement 777 if (cell->get_fe().n_dofs_per_face() == 0) 778 continue; 779 780 Assert(cell->face(face)->refinement_case() == 781 RefinementCase<dim - 1>::isotropic_refinement, 782 ExcNotImplemented()); 783 784 // in any case, faces can have at most two active fe indices, 785 // but here the face can have only one (namely the same as that 786 // from the cell we're sitting on), and each of the children can 787 // have only one as well. check this 788 AssertDimension(cell->face(face)->n_active_fe_indices(), 1); 789 Assert(cell->face(face)->fe_index_is_active( 790 cell->active_fe_index()) == true, 791 ExcInternalError()); 792 for (unsigned int c = 0; c < cell->face(face)->n_children(); 793 ++c) 794 if (!cell->neighbor_child_on_subface(face, c) 795 ->is_artificial()) 796 AssertDimension( 797 cell->face(face)->child(c)->n_active_fe_indices(), 1); 798 799 // right now, all that is implemented is the case that both 800 // sides use the same fe, and not only that but also that all 801 // lines bounding this face and the children have the same fe 802 for (unsigned int c = 0; c < cell->face(face)->n_children(); 803 ++c) 804 if (!cell->neighbor_child_on_subface(face, c) 805 ->is_artificial()) 806 { 807 Assert(cell->face(face)->child(c)->fe_index_is_active( 808 cell->active_fe_index()) == true, 809 ExcNotImplemented()); 810 for (unsigned int e = 0; e < 4; ++e) 811 { 812 Assert(cell->face(face) 813 ->child(c) 814 ->line(e) 815 ->n_active_fe_indices() == 1, 816 ExcNotImplemented()); 817 Assert(cell->face(face) 818 ->child(c) 819 ->line(e) 820 ->fe_index_is_active( 821 cell->active_fe_index()) == true, 822 ExcNotImplemented()); 823 } 824 } 825 for (unsigned int e = 0; e < 4; ++e) 826 { 827 Assert(cell->face(face)->line(e)->n_active_fe_indices() == 828 1, 829 ExcNotImplemented()); 830 Assert(cell->face(face)->line(e)->fe_index_is_active( 831 cell->active_fe_index()) == true, 832 ExcNotImplemented()); 833 } 834 835 // ok, start up the work 836 const FiniteElement<dim> &fe = cell->get_fe(); 837 const unsigned int fe_index = cell->active_fe_index(); 838 839 const unsigned int n_dofs_on_mother = fe.n_dofs_per_face(); 840 const unsigned int n_dofs_on_children = 841 (5 * fe.n_dofs_per_vertex() + 12 * fe.n_dofs_per_line() + 842 4 * fe.n_dofs_per_quad()); 843 844 // TODO[TL]: think about this and the following in case of 845 // anisotropic refinement 846 847 dofs_on_mother.resize(n_dofs_on_mother); 848 // we might not use all of those in case of artificial cells, so 849 // do not resize(), but reserve() and use push_back later. 850 dofs_on_children.clear(); 851 dofs_on_children.reserve(n_dofs_on_children); 852 853 Assert(n_dofs_on_mother == fe.constraints().n(), 854 ExcDimensionMismatch(n_dofs_on_mother, 855 fe.constraints().n())); 856 Assert(n_dofs_on_children == fe.constraints().m(), 857 ExcDimensionMismatch(n_dofs_on_children, 858 fe.constraints().m())); 859 860 const typename DoFHandler<dim_, spacedim>::face_iterator 861 this_face = cell->face(face); 862 863 // fill the dofs indices. Use same enumeration scheme as in 864 // @p{FiniteElement::constraints()} 865 unsigned int next_index = 0; 866 for (unsigned int vertex = 0; vertex < 4; ++vertex) 867 for (unsigned int dof = 0; dof != fe.n_dofs_per_vertex(); 868 ++dof) 869 dofs_on_mother[next_index++] = 870 this_face->vertex_dof_index(vertex, dof, fe_index); 871 for (unsigned int line = 0; line < 4; ++line) 872 for (unsigned int dof = 0; dof != fe.n_dofs_per_line(); ++dof) 873 dofs_on_mother[next_index++] = 874 this_face->line(line)->dof_index(dof, fe_index); 875 for (unsigned int dof = 0; dof != fe.n_dofs_per_quad(); ++dof) 876 dofs_on_mother[next_index++] = 877 this_face->dof_index(dof, fe_index); 878 AssertDimension(next_index, dofs_on_mother.size()); 879 880 // TODO: assert some consistency assumptions 881 882 // TODO[TL]: think about this in case of anisotropic refinement 883 884 Assert(dof_handler.get_triangulation() 885 .get_anisotropic_refinement_flag() || 886 ((this_face->child(0)->vertex_index(3) == 887 this_face->child(1)->vertex_index(2)) && 888 (this_face->child(0)->vertex_index(3) == 889 this_face->child(2)->vertex_index(1)) && 890 (this_face->child(0)->vertex_index(3) == 891 this_face->child(3)->vertex_index(0))), 892 ExcInternalError()); 893 894 for (unsigned int dof = 0; dof != fe.n_dofs_per_vertex(); ++dof) 895 dofs_on_children.push_back( 896 this_face->child(0)->vertex_dof_index(3, dof)); 897 898 // dof numbers on the centers of the lines bounding this face 899 for (unsigned int line = 0; line < 4; ++line) 900 for (unsigned int dof = 0; dof != fe.n_dofs_per_vertex(); 901 ++dof) 902 dofs_on_children.push_back( 903 this_face->line(line)->child(0)->vertex_dof_index( 904 1, dof, fe_index)); 905 906 // next the dofs on the lines interior to the face; the order of 907 // these lines is laid down in the FiniteElement class 908 // documentation 909 for (unsigned int dof = 0; dof < fe.n_dofs_per_line(); ++dof) 910 dofs_on_children.push_back( 911 this_face->child(0)->line(1)->dof_index(dof, fe_index)); 912 for (unsigned int dof = 0; dof < fe.n_dofs_per_line(); ++dof) 913 dofs_on_children.push_back( 914 this_face->child(2)->line(1)->dof_index(dof, fe_index)); 915 for (unsigned int dof = 0; dof < fe.n_dofs_per_line(); ++dof) 916 dofs_on_children.push_back( 917 this_face->child(0)->line(3)->dof_index(dof, fe_index)); 918 for (unsigned int dof = 0; dof < fe.n_dofs_per_line(); ++dof) 919 dofs_on_children.push_back( 920 this_face->child(1)->line(3)->dof_index(dof, fe_index)); 921 922 // dofs on the bordering lines 923 for (unsigned int line = 0; line < 4; ++line) 924 for (unsigned int child = 0; child < 2; ++child) 925 { 926 for (unsigned int dof = 0; dof != fe.n_dofs_per_line(); 927 ++dof) 928 dofs_on_children.push_back( 929 this_face->line(line)->child(child)->dof_index( 930 dof, fe_index)); 931 } 932 933 // finally, for the dofs interior to the four child faces 934 for (unsigned int child = 0; child < 4; ++child) 935 { 936 // skip artificial cells 937 if (cell->neighbor_child_on_subface(face, child) 938 ->is_artificial()) 939 continue; 940 for (unsigned int dof = 0; dof != fe.n_dofs_per_quad(); 941 ++dof) 942 dofs_on_children.push_back( 943 this_face->child(child)->dof_index(dof, fe_index)); 944 } 945 946 // note: can get fewer DoFs when we have artificial cells: 947 Assert(dofs_on_children.size() <= n_dofs_on_children, 948 ExcInternalError()); 949 950 // for each row in the AffineConstraints object for this line: 951 for (unsigned int row = 0; row != dofs_on_children.size(); 952 ++row) 953 { 954 constraints.add_line(dofs_on_children[row]); 955 for (unsigned int i = 0; i != dofs_on_mother.size(); ++i) 956 constraints.add_entry(dofs_on_children[row], 957 dofs_on_mother[i], 958 fe.constraints()(row, i)); 959 960 constraints.set_inhomogeneity(dofs_on_children[row], 0.); 961 } 962 } 963 else 964 { 965 // this face has no children, but it could still be that it is 966 // shared by two cells that use a different fe index. check a 967 // couple of things, but ignore the case that the neighbor is an 968 // artificial cell 969 if (!cell->at_boundary(face) && 970 !cell->neighbor(face)->is_artificial()) 971 { 972 Assert(cell->face(face)->n_active_fe_indices() == 1, 973 ExcNotImplemented()); 974 Assert(cell->face(face)->fe_index_is_active( 975 cell->active_fe_index()) == true, 976 ExcInternalError()); 977 } 978 } 979 } 980 } 981 982 983 984 template <int dim, int spacedim, typename number> 985 void make_hp_hanging_node_constraints(const DoFHandler<dim,spacedim> & dof_handler,AffineConstraints<number> & constraints)986 make_hp_hanging_node_constraints( 987 const DoFHandler<dim, spacedim> &dof_handler, 988 AffineConstraints<number> & constraints) 989 { 990 // note: this function is going to be hard to understand if you haven't 991 // read the hp paper. however, we try to follow the notation laid out 992 // there, so go read the paper before you try to understand what is going 993 // on here 994 995 996 // a matrix to be used for constraints below. declared here and simply 997 // resized down below to avoid permanent re-allocation of memory 998 FullMatrix<double> constraint_matrix; 999 1000 // similarly have arrays that will hold primary and dependent dof numbers, 1001 // as well as a scratch array needed for the complicated case below 1002 std::vector<types::global_dof_index> primary_dofs; 1003 std::vector<types::global_dof_index> dependent_dofs; 1004 std::vector<types::global_dof_index> scratch_dofs; 1005 1006 // caches for the face and subface interpolation matrices between 1007 // different (or the same) finite elements. we compute them only once, 1008 // namely the first time they are needed, and then just reuse them 1009 Table<2, std::unique_ptr<FullMatrix<double>>> face_interpolation_matrices( 1010 n_finite_elements(dof_handler), n_finite_elements(dof_handler)); 1011 Table<3, std::unique_ptr<FullMatrix<double>>> 1012 subface_interpolation_matrices( 1013 n_finite_elements(dof_handler), 1014 n_finite_elements(dof_handler), 1015 GeometryInfo<dim>::max_children_per_face); 1016 1017 // similarly have a cache for the matrices that are split into their 1018 // primary and dependent parts, and for which the primary part is 1019 // inverted. these two matrices are derived from the face interpolation 1020 // matrix 1021 // as described in the @ref hp_paper "hp paper" 1022 Table<2, 1023 std::unique_ptr<std::pair<FullMatrix<double>, FullMatrix<double>>>> 1024 split_face_interpolation_matrices(n_finite_elements(dof_handler), 1025 n_finite_elements(dof_handler)); 1026 1027 // finally, for each pair of finite elements, have a mask that states 1028 // which of the degrees of freedom on the coarse side of a refined face 1029 // will act as primary dofs. 1030 Table<2, std::unique_ptr<std::vector<bool>>> primary_dof_masks( 1031 n_finite_elements(dof_handler), n_finite_elements(dof_handler)); 1032 1033 // loop over all faces 1034 // 1035 // note that even though we may visit a face twice if the neighboring 1036 // cells are equally refined, we can only visit each face with hanging 1037 // nodes once 1038 for (const auto &cell : dof_handler.active_cell_iterators()) 1039 { 1040 // artificial cells can at best neighbor ghost cells, but we're not 1041 // interested in these interfaces 1042 if (cell->is_artificial()) 1043 continue; 1044 1045 for (const unsigned int face : cell->face_indices()) 1046 if (cell->face(face)->has_children()) 1047 { 1048 // first of all, make sure that we treat a case which is 1049 // possible, i.e. either no dofs on the face at all or no 1050 // anisotropic refinement 1051 if (cell->get_fe().n_dofs_per_face() == 0) 1052 continue; 1053 1054 Assert(cell->face(face)->refinement_case() == 1055 RefinementCase<dim - 1>::isotropic_refinement, 1056 ExcNotImplemented()); 1057 1058 // so now we've found a face of an active cell that has 1059 // children. that means that there are hanging nodes here. 1060 1061 // in any case, faces can have at most two sets of active fe 1062 // indices, but here the face can have only one (namely the same 1063 // as that from the cell we're sitting on), and each of the 1064 // children can have only one as well. check this 1065 Assert(cell->face(face)->n_active_fe_indices() == 1, 1066 ExcInternalError()); 1067 Assert(cell->face(face)->fe_index_is_active( 1068 cell->active_fe_index()) == true, 1069 ExcInternalError()); 1070 for (unsigned int c = 0; c < cell->face(face)->n_children(); 1071 ++c) 1072 if (!cell->neighbor_child_on_subface(face, c) 1073 ->is_artificial()) 1074 Assert(cell->face(face)->child(c)->n_active_fe_indices() == 1075 1, 1076 ExcInternalError()); 1077 1078 // first find out whether we can constrain each of the subfaces 1079 // to the mother face. in the lingo of the hp paper, this would 1080 // be the simple case. note that we can short-circuit this 1081 // decision if the dof_handler doesn't support hp at all 1082 // 1083 // ignore all interfaces with artificial cells 1084 FiniteElementDomination::Domination mother_face_dominates = 1085 FiniteElementDomination::either_element_can_dominate; 1086 1087 // auxiliary variable which holds FE indices of the mother face 1088 // and its subfaces. This knowledge will be needed in hp-case 1089 // with neither_element_dominates. 1090 std::set<unsigned int> fe_ind_face_subface; 1091 fe_ind_face_subface.insert(cell->active_fe_index()); 1092 1093 if (dof_handler.hp_capability_enabled) 1094 for (unsigned int c = 0; 1095 c < cell->face(face)->number_of_children(); 1096 ++c) 1097 { 1098 const auto subcell = 1099 cell->neighbor_child_on_subface(face, c); 1100 if (!subcell->is_artificial()) 1101 { 1102 mother_face_dominates = 1103 mother_face_dominates & 1104 (cell->get_fe().compare_for_domination( 1105 subcell->get_fe(), /*codim=*/1)); 1106 fe_ind_face_subface.insert( 1107 subcell->active_fe_index()); 1108 } 1109 } 1110 1111 switch (mother_face_dominates) 1112 { 1113 case FiniteElementDomination::this_element_dominates: 1114 case FiniteElementDomination::either_element_can_dominate: 1115 { 1116 // Case 1 (the simple case and the only case that can 1117 // happen for non-hp DoFHandlers): The coarse element 1118 // dominates the elements on the subfaces (or they are 1119 // all the same) 1120 // 1121 // so we are going to constrain the DoFs on the face 1122 // children against the DoFs on the face itself 1123 primary_dofs.resize(cell->get_fe().n_dofs_per_face()); 1124 1125 cell->face(face)->get_dof_indices( 1126 primary_dofs, cell->active_fe_index()); 1127 1128 // Now create constraints for the subfaces and 1129 // assemble it. ignore all interfaces with artificial 1130 // cells because we can only get to such interfaces if 1131 // the current cell is a ghost cell 1132 for (unsigned int c = 0; 1133 c < cell->face(face)->n_children(); 1134 ++c) 1135 { 1136 if (cell->neighbor_child_on_subface(face, c) 1137 ->is_artificial()) 1138 continue; 1139 1140 const typename DoFHandler<dim, spacedim>:: 1141 active_face_iterator subface = 1142 cell->face(face)->child(c); 1143 1144 Assert(subface->n_active_fe_indices() == 1, 1145 ExcInternalError()); 1146 1147 const unsigned int subface_fe_index = 1148 subface->nth_active_fe_index(0); 1149 1150 // we sometime run into the situation where for 1151 // example on one big cell we have a FE_Q(1) and on 1152 // the subfaces we have a mixture of FE_Q(1) and 1153 // FE_Nothing. In that case, the face domination is 1154 // either_element_can_dominate for the whole 1155 // collection of subfaces, but on the particular 1156 // subface between FE_Q(1) and FE_Nothing, there are 1157 // no constraints that we need to take care of. in 1158 // that case, just continue 1159 if (cell->get_fe().compare_for_domination( 1160 subface->get_fe(subface_fe_index), 1161 /*codim=*/1) == 1162 FiniteElementDomination::no_requirements) 1163 continue; 1164 1165 // Same procedure as for the mother cell. Extract 1166 // the face DoFs from the cell DoFs. 1167 dependent_dofs.resize( 1168 subface->get_fe(subface_fe_index) 1169 .n_dofs_per_face()); 1170 subface->get_dof_indices(dependent_dofs, 1171 subface_fe_index); 1172 1173 for (const types::global_dof_index dependent_dof : 1174 dependent_dofs) 1175 { 1176 (void)dependent_dof; 1177 Assert(dependent_dof != 1178 numbers::invalid_dof_index, 1179 ExcInternalError()); 1180 } 1181 1182 // Now create the element constraint for this 1183 // subface. 1184 // 1185 // As a side remark, one may wonder the following: 1186 // neighbor_child is clearly computed correctly, 1187 // i.e. taking into account face_orientation (just 1188 // look at the implementation of that function). 1189 // however, we don't care about this here, when we 1190 // ask for subface_interpolation on subface c. the 1191 // question rather is: do we have to translate 'c' 1192 // here as well? 1193 // 1194 // the answer is in fact 'no'. if one does that, 1195 // results are wrong: constraints are added twice 1196 // for the same pair of nodes but with differing 1197 // weights. in addition, one can look at the 1198 // deal.II/project_*_03 tests that look at exactly 1199 // this case: there, we have a mesh with at least 1200 // one face_orientation==false and hanging nodes, 1201 // and the results of those tests show that the 1202 // result of projection verifies the approximation 1203 // properties of a finite element onto that mesh 1204 ensure_existence_of_subface_matrix( 1205 cell->get_fe(), 1206 subface->get_fe(subface_fe_index), 1207 c, 1208 subface_interpolation_matrices 1209 [cell->active_fe_index()][subface_fe_index][c]); 1210 1211 // Add constraints to global AffineConstraints 1212 // object. 1213 filter_constraints(primary_dofs, 1214 dependent_dofs, 1215 *(subface_interpolation_matrices 1216 [cell->active_fe_index()] 1217 [subface_fe_index][c]), 1218 constraints); 1219 } // loop over subfaces 1220 1221 break; 1222 } // Case 1 1223 1224 case FiniteElementDomination::other_element_dominates: 1225 case FiniteElementDomination::neither_element_dominates: 1226 { 1227 // Case 2 (the "complex" case): at least one (the 1228 // neither_... case) of the finer elements or all of 1229 // them (the other_... case) is dominating. See the hp 1230 // paper for a way how to deal with this situation 1231 // 1232 // since this is something that can only happen for hp 1233 // dof handlers, add a check here... 1234 Assert(dof_handler.hp_capability_enabled == true, 1235 ExcInternalError()); 1236 1237 const dealii::hp::FECollection<dim, spacedim> 1238 &fe_collection = dof_handler.get_fe_collection(); 1239 // we first have to find the finite element that is able 1240 // to generate a space that all the other ones can be 1241 // constrained to. At this point we potentially have 1242 // different scenarios: 1243 // 1244 // 1) sub-faces dominate mother face and there is a 1245 // dominating FE among sub faces. We could loop over sub 1246 // faces to find the needed FE index. However, this will 1247 // not work in the case when ... 1248 // 1249 // 2) there is no dominating FE among sub faces (e.g. 1250 // Q1xQ2 vs Q2xQ1), but subfaces still dominate mother 1251 // face (e.g. Q2xQ2). To cover this case we would have 1252 // to find the least dominating element amongst all 1253 // finite elements on sub faces. 1254 // 1255 // 3) Finally, it could happen that we got here because 1256 // neither_element_dominates (e.g. Q1xQ1xQ2 and Q1xQ2xQ1 1257 // for subfaces and Q2xQ1xQ1 for mother face). This 1258 // requires finding the least dominating element amongst 1259 // all finite elements on sub faces and the mother face. 1260 // 1261 // Note that the last solution covers the first two 1262 // scenarios, thus we stick with it assuming that we 1263 // won't lose much time/efficiency. 1264 const unsigned int dominating_fe_index = 1265 fe_collection.find_dominating_fe_extended( 1266 fe_ind_face_subface, /*codim=*/1); 1267 1268 AssertThrow( 1269 dominating_fe_index != numbers::invalid_unsigned_int, 1270 ExcMessage( 1271 "Could not find a least face dominating FE.")); 1272 1273 const FiniteElement<dim, spacedim> &dominating_fe = 1274 dof_handler.get_fe(dominating_fe_index); 1275 1276 // first get the interpolation matrix from the mother to 1277 // the virtual dofs 1278 Assert(dominating_fe.n_dofs_per_face() <= 1279 cell->get_fe().n_dofs_per_face(), 1280 ExcInternalError()); 1281 1282 ensure_existence_of_face_matrix( 1283 dominating_fe, 1284 cell->get_fe(), 1285 face_interpolation_matrices[dominating_fe_index] 1286 [cell->active_fe_index()]); 1287 1288 // split this matrix into primary and dependent 1289 // components. invert the primary component 1290 ensure_existence_of_primary_dof_mask( 1291 cell->get_fe(), 1292 dominating_fe, 1293 (*face_interpolation_matrices 1294 [dominating_fe_index][cell->active_fe_index()]), 1295 primary_dof_masks[dominating_fe_index] 1296 [cell->active_fe_index()]); 1297 1298 ensure_existence_of_split_face_matrix( 1299 *face_interpolation_matrices[dominating_fe_index] 1300 [cell->active_fe_index()], 1301 (*primary_dof_masks[dominating_fe_index] 1302 [cell->active_fe_index()]), 1303 split_face_interpolation_matrices 1304 [dominating_fe_index][cell->active_fe_index()]); 1305 1306 const FullMatrix<double> 1307 &restrict_mother_to_virtual_primary_inv = 1308 (split_face_interpolation_matrices 1309 [dominating_fe_index][cell->active_fe_index()] 1310 ->first); 1311 1312 const FullMatrix<double> 1313 &restrict_mother_to_virtual_dependent = 1314 (split_face_interpolation_matrices 1315 [dominating_fe_index][cell->active_fe_index()] 1316 ->second); 1317 1318 // now compute the constraint matrix as the product 1319 // between the inverse matrix and the dependent part 1320 constraint_matrix.reinit( 1321 cell->get_fe().n_dofs_per_face() - 1322 dominating_fe.n_dofs_per_face(), 1323 dominating_fe.n_dofs_per_face()); 1324 restrict_mother_to_virtual_dependent.mmult( 1325 constraint_matrix, 1326 restrict_mother_to_virtual_primary_inv); 1327 1328 // then figure out the global numbers of primary and 1329 // dependent dofs and apply constraints 1330 scratch_dofs.resize(cell->get_fe().n_dofs_per_face()); 1331 cell->face(face)->get_dof_indices( 1332 scratch_dofs, cell->active_fe_index()); 1333 1334 // split dofs into primary and dependent components 1335 primary_dofs.clear(); 1336 dependent_dofs.clear(); 1337 for (unsigned int i = 0; 1338 i < cell->get_fe().n_dofs_per_face(); 1339 ++i) 1340 if ((*primary_dof_masks[dominating_fe_index] 1341 [cell 1342 ->active_fe_index()])[i] == 1343 true) 1344 primary_dofs.push_back(scratch_dofs[i]); 1345 else 1346 dependent_dofs.push_back(scratch_dofs[i]); 1347 1348 AssertDimension(primary_dofs.size(), 1349 dominating_fe.n_dofs_per_face()); 1350 AssertDimension(dependent_dofs.size(), 1351 cell->get_fe().n_dofs_per_face() - 1352 dominating_fe.n_dofs_per_face()); 1353 1354 filter_constraints(primary_dofs, 1355 dependent_dofs, 1356 constraint_matrix, 1357 constraints); 1358 1359 1360 1361 // next we have to deal with the subfaces. do as 1362 // discussed in the hp paper 1363 for (unsigned int sf = 0; 1364 sf < cell->face(face)->n_children(); 1365 ++sf) 1366 { 1367 // ignore interfaces with artificial cells as well 1368 // as interfaces between ghost cells in 2d 1369 if (cell->neighbor_child_on_subface(face, sf) 1370 ->is_artificial() || 1371 (dim == 2 && cell->is_ghost() && 1372 cell->neighbor_child_on_subface(face, sf) 1373 ->is_ghost())) 1374 continue; 1375 1376 Assert(cell->face(face) 1377 ->child(sf) 1378 ->n_active_fe_indices() == 1, 1379 ExcInternalError()); 1380 1381 const unsigned int subface_fe_index = 1382 cell->face(face)->child(sf)->nth_active_fe_index( 1383 0); 1384 const FiniteElement<dim, spacedim> &subface_fe = 1385 dof_handler.get_fe(subface_fe_index); 1386 1387 // first get the interpolation matrix from the 1388 // subface to the virtual dofs 1389 Assert(dominating_fe.n_dofs_per_face() <= 1390 subface_fe.n_dofs_per_face(), 1391 ExcInternalError()); 1392 ensure_existence_of_subface_matrix( 1393 dominating_fe, 1394 subface_fe, 1395 sf, 1396 subface_interpolation_matrices 1397 [dominating_fe_index][subface_fe_index][sf]); 1398 1399 const FullMatrix<double> 1400 &restrict_subface_to_virtual = *( 1401 subface_interpolation_matrices 1402 [dominating_fe_index][subface_fe_index][sf]); 1403 1404 constraint_matrix.reinit( 1405 subface_fe.n_dofs_per_face(), 1406 dominating_fe.n_dofs_per_face()); 1407 1408 restrict_subface_to_virtual.mmult( 1409 constraint_matrix, 1410 restrict_mother_to_virtual_primary_inv); 1411 1412 dependent_dofs.resize(subface_fe.n_dofs_per_face()); 1413 cell->face(face)->child(sf)->get_dof_indices( 1414 dependent_dofs, subface_fe_index); 1415 1416 filter_constraints(primary_dofs, 1417 dependent_dofs, 1418 constraint_matrix, 1419 constraints); 1420 } // loop over subfaces 1421 1422 break; 1423 } // Case 2 1424 1425 case FiniteElementDomination::no_requirements: 1426 // there are no continuity requirements between the two 1427 // elements. record no constraints 1428 break; 1429 1430 default: 1431 // we shouldn't get here 1432 Assert(false, ExcInternalError()); 1433 } 1434 } 1435 else 1436 { 1437 // this face has no children, but it could still be that it is 1438 // shared by two cells that use a different fe index 1439 Assert(cell->face(face)->fe_index_is_active( 1440 cell->active_fe_index()) == true, 1441 ExcInternalError()); 1442 1443 // see if there is a neighbor that is an artificial cell. in 1444 // that case, we're not interested in this interface. we test 1445 // this case first since artificial cells may not have an 1446 // active_fe_index set, etc 1447 if (!cell->at_boundary(face) && 1448 cell->neighbor(face)->is_artificial()) 1449 continue; 1450 1451 // Only if there is a neighbor with a different active_fe_index 1452 // and the same h-level, some action has to be taken. 1453 if ((dof_handler.hp_capability_enabled) && 1454 !cell->face(face)->at_boundary() && 1455 (cell->neighbor(face)->active_fe_index() != 1456 cell->active_fe_index()) && 1457 (!cell->face(face)->has_children() && 1458 !cell->neighbor_is_coarser(face))) 1459 { 1460 const typename DoFHandler<dim, 1461 spacedim>::level_cell_iterator 1462 neighbor = cell->neighbor(face); 1463 1464 // see which side of the face we have to constrain 1465 switch ( 1466 cell->get_fe().compare_for_domination(neighbor->get_fe(), 1467 /*codim=*/1)) 1468 { 1469 case FiniteElementDomination::this_element_dominates: 1470 { 1471 // Get DoFs on dominating and dominated side of the 1472 // face 1473 primary_dofs.resize( 1474 cell->get_fe().n_dofs_per_face()); 1475 cell->face(face)->get_dof_indices( 1476 primary_dofs, cell->active_fe_index()); 1477 1478 // break if the n_primary_dofs == 0, because we are 1479 // attempting to constrain to an element that has no 1480 // face dofs 1481 if (primary_dofs.size() == 0) 1482 break; 1483 1484 dependent_dofs.resize( 1485 neighbor->get_fe().n_dofs_per_face()); 1486 cell->face(face)->get_dof_indices( 1487 dependent_dofs, neighbor->active_fe_index()); 1488 1489 // make sure the element constraints for this face 1490 // are available 1491 ensure_existence_of_face_matrix( 1492 cell->get_fe(), 1493 neighbor->get_fe(), 1494 face_interpolation_matrices 1495 [cell->active_fe_index()] 1496 [neighbor->active_fe_index()]); 1497 1498 // Add constraints to global constraint matrix. 1499 filter_constraints( 1500 primary_dofs, 1501 dependent_dofs, 1502 *(face_interpolation_matrices 1503 [cell->active_fe_index()] 1504 [neighbor->active_fe_index()]), 1505 constraints); 1506 1507 break; 1508 } 1509 1510 case FiniteElementDomination::other_element_dominates: 1511 { 1512 // we don't do anything here since we will come back 1513 // to this face from the other cell, at which time 1514 // we will fall into the first case clause above 1515 break; 1516 } 1517 1518 case FiniteElementDomination:: 1519 either_element_can_dominate: 1520 { 1521 // it appears as if neither element has any 1522 // constraints on its neighbor. this may be because 1523 // neither element has any DoFs on faces at all. or 1524 // that the two elements are actually the same, 1525 // although they happen to run under different 1526 // fe_indices (this is what happens in 1527 // hp/hp_hanging_nodes_01 for example). 1528 // 1529 // another possibility is what happens in crash_13. 1530 // there, we have FESystem(FE_Q(1),FE_DGQ(0)) vs. 1531 // FESystem(FE_Q(1),FE_DGQ(1)). neither of them 1532 // dominates the other. 1533 // 1534 // a final possibility is that we have something 1535 // like FESystem(FE_Q(1),FE_Q(1)) vs 1536 // FESystem(FE_Q(1),FE_Nothing()), see 1537 // hp/fe_nothing_18/19. 1538 // 1539 // in any case, the point is that it doesn't matter. 1540 // there is nothing to do here. 1541 break; 1542 } 1543 1544 case FiniteElementDomination::neither_element_dominates: 1545 { 1546 // make sure we don't get here twice from each cell 1547 if (cell < neighbor) 1548 break; 1549 1550 // our best bet is to find the common space among 1551 // other FEs in FECollection and then constrain both 1552 // FEs to that one. More precisely, we follow the 1553 // strategy outlined on page 17 of the hp paper: 1554 // First we find the dominant FE space S. Then we 1555 // divide our dofs in primary and dependent such 1556 // that I^{face,primary}_{S^{face}->S} is 1557 // invertible. And finally constrain dependent dofs 1558 // to primary dofs based on the interpolation 1559 // matrix. 1560 1561 const unsigned int this_fe_index = 1562 cell->active_fe_index(); 1563 const unsigned int neighbor_fe_index = 1564 neighbor->active_fe_index(); 1565 std::set<unsigned int> fes; 1566 fes.insert(this_fe_index); 1567 fes.insert(neighbor_fe_index); 1568 const dealii::hp::FECollection<dim, spacedim> 1569 &fe_collection = dof_handler.get_fe_collection(); 1570 1571 const unsigned int dominating_fe_index = 1572 fe_collection.find_dominating_fe_extended( 1573 fes, /*codim=*/1); 1574 1575 AssertThrow( 1576 dominating_fe_index != 1577 numbers::invalid_unsigned_int, 1578 ExcMessage( 1579 "Could not find the dominating FE for " + 1580 cell->get_fe().get_name() + " and " + 1581 neighbor->get_fe().get_name() + 1582 " inside FECollection.")); 1583 1584 const FiniteElement<dim, spacedim> &dominating_fe = 1585 fe_collection[dominating_fe_index]; 1586 1587 // TODO: until we hit the second face, the code is a 1588 // copy-paste from h-refinement case... 1589 1590 // first get the interpolation matrix from main FE 1591 // to the virtual dofs 1592 Assert(dominating_fe.n_dofs_per_face() <= 1593 cell->get_fe().n_dofs_per_face(), 1594 ExcInternalError()); 1595 1596 ensure_existence_of_face_matrix( 1597 dominating_fe, 1598 cell->get_fe(), 1599 face_interpolation_matrices 1600 [dominating_fe_index][cell->active_fe_index()]); 1601 1602 // split this matrix into primary and dependent 1603 // components. invert the primary component 1604 ensure_existence_of_primary_dof_mask( 1605 cell->get_fe(), 1606 dominating_fe, 1607 (*face_interpolation_matrices 1608 [dominating_fe_index] 1609 [cell->active_fe_index()]), 1610 primary_dof_masks[dominating_fe_index] 1611 [cell->active_fe_index()]); 1612 1613 ensure_existence_of_split_face_matrix( 1614 *face_interpolation_matrices 1615 [dominating_fe_index][cell->active_fe_index()], 1616 (*primary_dof_masks[dominating_fe_index] 1617 [cell->active_fe_index()]), 1618 split_face_interpolation_matrices 1619 [dominating_fe_index][cell->active_fe_index()]); 1620 1621 const FullMatrix< 1622 double> &restrict_mother_to_virtual_primary_inv = 1623 (split_face_interpolation_matrices 1624 [dominating_fe_index][cell->active_fe_index()] 1625 ->first); 1626 1627 const FullMatrix< 1628 double> &restrict_mother_to_virtual_dependent = 1629 (split_face_interpolation_matrices 1630 [dominating_fe_index][cell->active_fe_index()] 1631 ->second); 1632 1633 // now compute the constraint matrix as the product 1634 // between the inverse matrix and the dependent part 1635 constraint_matrix.reinit( 1636 cell->get_fe().n_dofs_per_face() - 1637 dominating_fe.n_dofs_per_face(), 1638 dominating_fe.n_dofs_per_face()); 1639 restrict_mother_to_virtual_dependent.mmult( 1640 constraint_matrix, 1641 restrict_mother_to_virtual_primary_inv); 1642 1643 // then figure out the global numbers of primary and 1644 // dependent dofs and apply constraints 1645 scratch_dofs.resize( 1646 cell->get_fe().n_dofs_per_face()); 1647 cell->face(face)->get_dof_indices( 1648 scratch_dofs, cell->active_fe_index()); 1649 1650 // split dofs into primary and dependent components 1651 primary_dofs.clear(); 1652 dependent_dofs.clear(); 1653 for (unsigned int i = 0; 1654 i < cell->get_fe().n_dofs_per_face(); 1655 ++i) 1656 if ((*primary_dof_masks[dominating_fe_index] 1657 [cell->active_fe_index()]) 1658 [i] == true) 1659 primary_dofs.push_back(scratch_dofs[i]); 1660 else 1661 dependent_dofs.push_back(scratch_dofs[i]); 1662 1663 AssertDimension(primary_dofs.size(), 1664 dominating_fe.n_dofs_per_face()); 1665 AssertDimension(dependent_dofs.size(), 1666 cell->get_fe().n_dofs_per_face() - 1667 dominating_fe.n_dofs_per_face()); 1668 1669 filter_constraints(primary_dofs, 1670 dependent_dofs, 1671 constraint_matrix, 1672 constraints); 1673 1674 // now do the same for another FE this is pretty 1675 // much the same we do above to resolve h-refinement 1676 // constraints 1677 Assert(dominating_fe.n_dofs_per_face() <= 1678 neighbor->get_fe().n_dofs_per_face(), 1679 ExcInternalError()); 1680 1681 ensure_existence_of_face_matrix( 1682 dominating_fe, 1683 neighbor->get_fe(), 1684 face_interpolation_matrices 1685 [dominating_fe_index] 1686 [neighbor->active_fe_index()]); 1687 1688 const FullMatrix<double> 1689 &restrict_secondface_to_virtual = 1690 *(face_interpolation_matrices 1691 [dominating_fe_index] 1692 [neighbor->active_fe_index()]); 1693 1694 constraint_matrix.reinit( 1695 neighbor->get_fe().n_dofs_per_face(), 1696 dominating_fe.n_dofs_per_face()); 1697 1698 restrict_secondface_to_virtual.mmult( 1699 constraint_matrix, 1700 restrict_mother_to_virtual_primary_inv); 1701 1702 dependent_dofs.resize( 1703 neighbor->get_fe().n_dofs_per_face()); 1704 cell->face(face)->get_dof_indices( 1705 dependent_dofs, neighbor->active_fe_index()); 1706 1707 filter_constraints(primary_dofs, 1708 dependent_dofs, 1709 constraint_matrix, 1710 constraints); 1711 1712 break; 1713 } 1714 1715 case FiniteElementDomination::no_requirements: 1716 { 1717 // nothing to do here 1718 break; 1719 } 1720 1721 default: 1722 // we shouldn't get here 1723 Assert(false, ExcInternalError()); 1724 } 1725 } 1726 } 1727 } 1728 } 1729 } // namespace internal 1730 1731 1732 1733 template <int dim, int spacedim, typename number> 1734 void make_hanging_node_constraints(const DoFHandler<dim,spacedim> & dof_handler,AffineConstraints<number> & constraints)1735 make_hanging_node_constraints(const DoFHandler<dim, spacedim> &dof_handler, 1736 AffineConstraints<number> & constraints) 1737 { 1738 // Decide whether to use the new or old make_hanging_node_constraints 1739 // function. If all the FiniteElement or all elements in a FECollection 1740 // support the new face constraint matrix, the new code will be used. 1741 // Otherwise, the old implementation is used for the moment. 1742 if (dof_handler.get_fe_collection().hp_constraints_are_implemented()) 1743 internal::make_hp_hanging_node_constraints(dof_handler, constraints); 1744 else 1745 internal::make_oldstyle_hanging_node_constraints( 1746 dof_handler, constraints, std::integral_constant<int, dim>()); 1747 } 1748 1749 1750 1751 namespace 1752 { 1753 /** 1754 * @internal 1755 * 1756 * Internally used in make_periodicity_constraints. 1757 * 1758 * enter constraints for periodicity into the given AffineConstraints 1759 * object. this function is called when at least one of the two face 1760 * iterators corresponds to an active object without further children 1761 * 1762 * @param transformation A matrix that maps degrees of freedom from one face 1763 * to another. If the DoFs on the two faces are supposed to match exactly, 1764 * then the matrix so provided will be the identity matrix. if face 2 is 1765 * once refined from face 1, then the matrix needs to be the interpolation 1766 * matrix from a face to this particular child 1767 * 1768 * @precondition: face_1 is supposed to be active 1769 * 1770 * @note We have to be careful not to accidentally create constraint 1771 * cycles when adding periodic constraints: For example, as the 1772 * corresponding testcase bits/periodicity_05 demonstrates, we can 1773 * occasionally get into trouble if we already have the constraint 1774 * x1 == x2 and want to insert x2 == x1. We avoid this by skipping 1775 * such "identity constraints" if the opposite constraint already 1776 * exists. 1777 */ 1778 template <typename FaceIterator, typename number> 1779 void set_periodicity_constraints(const FaceIterator & face_1,const typename identity<FaceIterator>::type & face_2,const FullMatrix<double> & transformation,AffineConstraints<number> & affine_constraints,const ComponentMask & component_mask,const bool face_orientation,const bool face_flip,const bool face_rotation,const number periodicity_factor)1780 set_periodicity_constraints( 1781 const FaceIterator & face_1, 1782 const typename identity<FaceIterator>::type &face_2, 1783 const FullMatrix<double> & transformation, 1784 AffineConstraints<number> & affine_constraints, 1785 const ComponentMask & component_mask, 1786 const bool face_orientation, 1787 const bool face_flip, 1788 const bool face_rotation, 1789 const number periodicity_factor) 1790 { 1791 static const int dim = FaceIterator::AccessorType::dimension; 1792 static const int spacedim = FaceIterator::AccessorType::space_dimension; 1793 1794 // we should be in the case where face_1 is active, i.e. has no children: 1795 Assert(!face_1->has_children(), ExcInternalError()); 1796 1797 Assert(face_1->n_active_fe_indices() == 1, ExcInternalError()); 1798 1799 // If face_2 does have children, then we need to iterate over these 1800 // children and set periodic constraints in the inverse direction: 1801 1802 if (face_2->has_children()) 1803 { 1804 Assert(face_2->n_children() == 1805 GeometryInfo<dim>::max_children_per_face, 1806 ExcNotImplemented()); 1807 1808 const unsigned int dofs_per_face = 1809 face_1->get_fe(face_1->nth_active_fe_index(0)).n_dofs_per_face(); 1810 FullMatrix<double> child_transformation(dofs_per_face, dofs_per_face); 1811 FullMatrix<double> subface_interpolation(dofs_per_face, 1812 dofs_per_face); 1813 1814 for (unsigned int c = 0; c < face_2->n_children(); ++c) 1815 { 1816 // get the interpolation matrix recursively from the one that 1817 // interpolated from face_1 to face_2 by multiplying from the left 1818 // with the one that interpolates from face_2 to its child 1819 const auto &fe = face_1->get_fe(face_1->nth_active_fe_index(0)); 1820 fe.get_subface_interpolation_matrix(fe, c, subface_interpolation); 1821 subface_interpolation.mmult(child_transformation, transformation); 1822 1823 set_periodicity_constraints(face_1, 1824 face_2->child(c), 1825 child_transformation, 1826 affine_constraints, 1827 component_mask, 1828 face_orientation, 1829 face_flip, 1830 face_rotation, 1831 periodicity_factor); 1832 } 1833 return; 1834 } 1835 1836 // 1837 // If we reached this point then both faces are active. Now all 1838 // that is left is to match the corresponding DoFs of both faces. 1839 // 1840 1841 const unsigned int face_1_index = face_1->nth_active_fe_index(0); 1842 const unsigned int face_2_index = face_2->nth_active_fe_index(0); 1843 Assert(face_1->get_fe(face_1_index) == face_2->get_fe(face_2_index), 1844 ExcMessage( 1845 "Matching periodic cells need to use the same finite element")); 1846 1847 const FiniteElement<dim, spacedim> &fe = face_1->get_fe(face_1_index); 1848 1849 Assert(component_mask.represents_n_components(fe.n_components()), 1850 ExcMessage( 1851 "The number of components in the mask has to be either " 1852 "zero or equal to the number of components in the finite " 1853 "element.")); 1854 1855 const unsigned int dofs_per_face = fe.n_dofs_per_face(); 1856 1857 std::vector<types::global_dof_index> dofs_1(dofs_per_face); 1858 std::vector<types::global_dof_index> dofs_2(dofs_per_face); 1859 1860 face_1->get_dof_indices(dofs_1, face_1_index); 1861 face_2->get_dof_indices(dofs_2, face_2_index); 1862 1863 // If either of the two faces has an invalid dof index, stop. This is 1864 // so that there is no attempt to match artificial cells of parallel 1865 // distributed triangulations. 1866 // 1867 // While it seems like we ought to be able to avoid even calling 1868 // set_periodicity_constraints for artificial faces, this situation 1869 // can arise when a face that is being made periodic is only 1870 // partially touched by the local subdomain. 1871 // make_periodicity_constraints will be called recursively even for 1872 // the section of the face that is not touched by the local 1873 // subdomain. 1874 // 1875 // Until there is a better way to determine if the cells that 1876 // neighbor a face are artificial, we simply test to see if the face 1877 // does not have a valid dof initialization. 1878 1879 for (unsigned int i = 0; i < dofs_per_face; i++) 1880 if (dofs_1[i] == numbers::invalid_dof_index || 1881 dofs_2[i] == numbers::invalid_dof_index) 1882 { 1883 return; 1884 } 1885 1886 // Well, this is a hack: 1887 // 1888 // There is no 1889 // face_to_face_index(face_index, 1890 // face_orientation, 1891 // face_flip, 1892 // face_rotation) 1893 // function in FiniteElementData, so we have to use 1894 // face_to_cell_index(face_index, face 1895 // face_orientation, 1896 // face_flip, 1897 // face_rotation) 1898 // But this will give us an index on a cell - something we cannot work 1899 // with directly. But luckily we can match them back :-] 1900 1901 std::map<unsigned int, unsigned int> cell_to_rotated_face_index; 1902 1903 // Build up a cell to face index for face_2: 1904 for (unsigned int i = 0; i < dofs_per_face; ++i) 1905 { 1906 const unsigned int cell_index = 1907 fe.face_to_cell_index(i, 1908 0, /* It doesn't really matter, just 1909 * assume we're on the first face... 1910 */ 1911 true, 1912 false, 1913 false // default orientation 1914 ); 1915 cell_to_rotated_face_index[cell_index] = i; 1916 } 1917 1918 // 1919 // Loop over all dofs on face 2 and constrain them against all 1920 // matching dofs on face 1: 1921 // 1922 1923 for (unsigned int i = 0; i < dofs_per_face; ++i) 1924 { 1925 // Obey the component mask 1926 if ((component_mask.n_selected_components(fe.n_components()) != 1927 fe.n_components()) && 1928 !component_mask[fe.face_system_to_component_index(i).first]) 1929 continue; 1930 1931 // We have to be careful to treat so called "identity 1932 // constraints" special. These are constraints of the form 1933 // x1 == constraint_factor * x_2. In this case, if the constraint 1934 // x2 == 1./constraint_factor * x1 already exists we are in trouble. 1935 // 1936 // Consequently, we have to check that we have indeed such an 1937 // "identity constraint". We do this by looping over all entries 1938 // of the row of the transformation matrix and check whether we 1939 // find exactly one nonzero entry. If this is the case, set 1940 // "is_identity_constrained" to true and record the corresponding 1941 // index and constraint_factor. 1942 1943 bool is_identity_constrained = false; 1944 unsigned int target = numbers::invalid_unsigned_int; 1945 number constraint_factor = periodicity_factor; 1946 1947 constexpr double eps = 1.e-13; 1948 for (unsigned int jj = 0; jj < dofs_per_face; ++jj) 1949 { 1950 const auto entry = transformation(i, jj); 1951 if (std::abs(entry) > eps) 1952 { 1953 if (is_identity_constrained) 1954 { 1955 // We did encounter more than one nonzero entry, so 1956 // the dof is not identity constrained. Set the 1957 // boolean to false and break out of the for loop. 1958 is_identity_constrained = false; 1959 break; 1960 } 1961 is_identity_constrained = true; 1962 target = jj; 1963 constraint_factor = entry * periodicity_factor; 1964 } 1965 } 1966 1967 // Next, we work on all constraints that are not identity 1968 // constraints, i.e., constraints that involve an interpolation 1969 // step that constrains the current dof (on face 2) to more than 1970 // one dof on face 1. 1971 1972 if (!is_identity_constrained) 1973 { 1974 // The current dof is already constrained. There is nothing 1975 // left to do. 1976 if (affine_constraints.is_constrained(dofs_2[i])) 1977 continue; 1978 1979 // Enter the constraint piece by piece: 1980 affine_constraints.add_line(dofs_2[i]); 1981 1982 for (unsigned int jj = 0; jj < dofs_per_face; ++jj) 1983 { 1984 // Get the correct dof index on face_1 respecting the 1985 // given orientation: 1986 const unsigned int j = 1987 cell_to_rotated_face_index[fe.face_to_cell_index( 1988 jj, 0, face_orientation, face_flip, face_rotation)]; 1989 1990 if (std::abs(transformation(i, jj)) > eps) 1991 affine_constraints.add_entry(dofs_2[i], 1992 dofs_1[j], 1993 transformation(i, jj)); 1994 } 1995 1996 // Continue with next dof. 1997 continue; 1998 } 1999 2000 // We are left with an "identity constraint". 2001 2002 // Get the correct dof index on face_1 respecting the given 2003 // orientation: 2004 const unsigned int j = 2005 cell_to_rotated_face_index[fe.face_to_cell_index( 2006 target, 0, face_orientation, face_flip, face_rotation)]; 2007 2008 auto dof_left = dofs_1[j]; 2009 auto dof_right = dofs_2[i]; 2010 2011 // If dof_left is already constrained, or dof_left < dof_right we 2012 // flip the order to ensure that dofs are constrained in a stable 2013 // manner on different MPI processes. 2014 if (affine_constraints.is_constrained(dof_left) || 2015 (dof_left < dof_right && 2016 !affine_constraints.is_constrained(dof_right))) 2017 { 2018 std::swap(dof_left, dof_right); 2019 constraint_factor = 1. / constraint_factor; 2020 } 2021 2022 // Next, we try to enter the constraint 2023 // dof_left = constraint_factor * dof_right; 2024 2025 // If both degrees of freedom are constrained, there is nothing we 2026 // can do. Simply continue with the next dof. 2027 if (affine_constraints.is_constrained(dof_left) && 2028 affine_constraints.is_constrained(dof_right)) 2029 continue; 2030 2031 // We have to be careful that adding the current identity 2032 // constraint does not create a constraint cycle. Thus, check for 2033 // a dependency cycle: 2034 2035 bool constraints_are_cyclic = true; 2036 number cycle_constraint_factor = constraint_factor; 2037 2038 for (auto test_dof = dof_right; test_dof != dof_left;) 2039 { 2040 if (!affine_constraints.is_constrained(test_dof)) 2041 { 2042 constraints_are_cyclic = false; 2043 break; 2044 } 2045 2046 const auto &constraint_entries = 2047 *affine_constraints.get_constraint_entries(test_dof); 2048 if (constraint_entries.size() == 1) 2049 { 2050 test_dof = constraint_entries[0].first; 2051 cycle_constraint_factor *= constraint_entries[0].second; 2052 } 2053 else 2054 { 2055 constraints_are_cyclic = false; 2056 break; 2057 } 2058 } 2059 2060 // In case of a dependency cycle we, either 2061 // - do nothing if cycle_constraint_factor == 1. In this case all 2062 // degrees 2063 // of freedom are already periodically constrained, 2064 // - otherwise, force all dofs to zero (by setting dof_left to 2065 // zero). The reasoning behind this is the fact that 2066 // cycle_constraint_factor != 1 occurs in situations such as 2067 // x1 == x2 and x2 == -1. * x1. This system is only solved by 2068 // x_1 = x_2 = 0. 2069 2070 if (constraints_are_cyclic) 2071 { 2072 if (std::abs(cycle_constraint_factor - 1.) > eps) 2073 affine_constraints.add_line(dof_left); 2074 } 2075 else 2076 { 2077 affine_constraints.add_line(dof_left); 2078 affine_constraints.add_entry(dof_left, 2079 dof_right, 2080 constraint_factor); 2081 // The number 1e10 in the assert below is arbitrary. If the 2082 // absolute value of constraint_factor is too large, then probably 2083 // the absolute value of periodicity_factor is too large or too 2084 // small. This would be equivalent to an evanescent wave that has 2085 // a very small wavelength. A quick calculation shows that if 2086 // |periodicity_factor| > 1e10 -> |np.exp(ikd)|> 1e10, therefore k 2087 // is imaginary (evanescent wave) and the evanescent wavelength is 2088 // 0.27 times smaller than the dimension of the structure, 2089 // lambda=((2*pi)/log(1e10))*d. Imaginary wavenumbers can be 2090 // interesting in some cases 2091 // (https://doi.org/10.1103/PhysRevA.94.033813).In order to 2092 // implement the case of in which the wavevector can be imaginary 2093 // it would be necessary to rewrite this function and the dof 2094 // ordering method should be modified. 2095 // Let's take the following constraint a*x1 + b*x2 = 0. You could 2096 // just always pick x1 = b/a*x2, but in practice this is not so 2097 // stable if a could be a small number -- intended to be zero, but 2098 // just very small due to roundoff. Of course, constraining x2 in 2099 // terms of x1 has the same problem. So one chooses x1 = b/a*x2 if 2100 // |b|<|a|, and x2 = a/b*x1 if |a|<|b|. 2101 Assert( 2102 std::abs(constraint_factor) < 1e10, 2103 ExcMessage( 2104 "The periodicity constraint is too large. The parameter periodicity_factor might be too large or too small.")); 2105 } 2106 } /* for dofs_per_face */ 2107 } 2108 2109 2110 2111 // Internally used in make_periodicity_constraints. 2112 // 2113 // Build up a (possibly rotated) interpolation matrix that is used in 2114 // set_periodicity_constraints with the help of user supplied matrix and 2115 // first_vector_components. 2116 template <int dim, int spacedim> 2117 FullMatrix<double> compute_transformation(const FiniteElement<dim,spacedim> & fe,const FullMatrix<double> & matrix,const std::vector<unsigned int> & first_vector_components)2118 compute_transformation( 2119 const FiniteElement<dim, spacedim> &fe, 2120 const FullMatrix<double> & matrix, 2121 const std::vector<unsigned int> & first_vector_components) 2122 { 2123 Assert(matrix.m() == matrix.n(), ExcInternalError()); 2124 2125 const unsigned int n_dofs_per_face = fe.n_dofs_per_face(); 2126 2127 if (matrix.m() == n_dofs_per_face) 2128 { 2129 // In case of m == n == n_dofs_per_face the supplied matrix is already 2130 // an interpolation matrix, so we use it directly: 2131 return matrix; 2132 } 2133 2134 if (first_vector_components.empty() && matrix.m() == 0) 2135 { 2136 // Just the identity matrix in case no rotation is specified: 2137 return IdentityMatrix(n_dofs_per_face); 2138 } 2139 2140 // The matrix describes a rotation and we have to build a transformation 2141 // matrix, we assume that for a 0* rotation we would have to build the 2142 // identity matrix 2143 2144 Assert(matrix.m() == spacedim, ExcInternalError()) 2145 2146 Quadrature<dim - 1> 2147 quadrature(fe.get_unit_face_support_points()); 2148 2149 // have an array that stores the location of each vector-dof tuple we want 2150 // to rotate. 2151 using DoFTuple = std::array<unsigned int, spacedim>; 2152 2153 // start with a pristine interpolation matrix... 2154 FullMatrix<double> transformation = IdentityMatrix(n_dofs_per_face); 2155 2156 for (unsigned int i = 0; i < n_dofs_per_face; ++i) 2157 { 2158 std::vector<unsigned int>::const_iterator comp_it = 2159 std::find(first_vector_components.begin(), 2160 first_vector_components.end(), 2161 fe.face_system_to_component_index(i).first); 2162 if (comp_it != first_vector_components.end()) 2163 { 2164 const unsigned int first_vector_component = *comp_it; 2165 2166 // find corresponding other components of vector 2167 DoFTuple vector_dofs; 2168 vector_dofs[0] = i; 2169 unsigned int n_found = 1; 2170 2171 Assert( 2172 *comp_it + spacedim <= fe.n_components(), 2173 ExcMessage( 2174 "Error: the finite element does not have enough components " 2175 "to define rotated periodic boundaries.")); 2176 2177 for (unsigned int k = 0; k < n_dofs_per_face; ++k) 2178 if ((k != i) && (quadrature.point(k) == quadrature.point(i)) && 2179 (fe.face_system_to_component_index(k).first >= 2180 first_vector_component) && 2181 (fe.face_system_to_component_index(k).first < 2182 first_vector_component + spacedim)) 2183 { 2184 vector_dofs[fe.face_system_to_component_index(k).first - 2185 first_vector_component] = k; 2186 n_found++; 2187 if (n_found == dim) 2188 break; 2189 } 2190 2191 // ... and rotate all dofs belonging to vector valued components 2192 // that are selected by first_vector_components: 2193 for (int i = 0; i < spacedim; ++i) 2194 { 2195 transformation[vector_dofs[i]][vector_dofs[i]] = 0.; 2196 for (int j = 0; j < spacedim; ++j) 2197 transformation[vector_dofs[i]][vector_dofs[j]] = 2198 matrix[i][j]; 2199 } 2200 } 2201 } 2202 return transformation; 2203 } 2204 } /*namespace*/ 2205 2206 2207 // Low level interface: 2208 2209 2210 template <typename FaceIterator, typename number> 2211 void make_periodicity_constraints(const FaceIterator & face_1,const typename identity<FaceIterator>::type & face_2,AffineConstraints<number> & affine_constraints,const ComponentMask & component_mask,const bool face_orientation,const bool face_flip,const bool face_rotation,const FullMatrix<double> & matrix,const std::vector<unsigned int> & first_vector_components,const number periodicity_factor)2212 make_periodicity_constraints( 2213 const FaceIterator & face_1, 2214 const typename identity<FaceIterator>::type &face_2, 2215 AffineConstraints<number> & affine_constraints, 2216 const ComponentMask & component_mask, 2217 const bool face_orientation, 2218 const bool face_flip, 2219 const bool face_rotation, 2220 const FullMatrix<double> & matrix, 2221 const std::vector<unsigned int> & first_vector_components, 2222 const number periodicity_factor) 2223 { 2224 static const int dim = FaceIterator::AccessorType::dimension; 2225 static const int spacedim = FaceIterator::AccessorType::space_dimension; 2226 2227 Assert((dim != 1) || (face_orientation == true && face_flip == false && 2228 face_rotation == false), 2229 ExcMessage("The supplied orientation " 2230 "(face_orientation, face_flip, face_rotation) " 2231 "is invalid for 1D")); 2232 2233 Assert((dim != 2) || (face_orientation == true && face_rotation == false), 2234 ExcMessage("The supplied orientation " 2235 "(face_orientation, face_flip, face_rotation) " 2236 "is invalid for 2D")); 2237 2238 Assert(face_1 != face_2, 2239 ExcMessage("face_1 and face_2 are equal! Cannot constrain DoFs " 2240 "on the very same face")); 2241 2242 Assert(face_1->at_boundary() && face_2->at_boundary(), 2243 ExcMessage("Faces for periodicity constraints must be on the " 2244 "boundary")); 2245 2246 Assert(matrix.m() == matrix.n(), 2247 ExcMessage("The supplied (rotation or interpolation) matrix must " 2248 "be a square matrix")); 2249 2250 Assert(first_vector_components.empty() || matrix.m() == spacedim, 2251 ExcMessage("first_vector_components is nonempty, so matrix must " 2252 "be a rotation matrix exactly of size spacedim")); 2253 2254 #ifdef DEBUG 2255 if (!face_1->has_children()) 2256 { 2257 Assert(face_1->n_active_fe_indices() == 1, ExcInternalError()); 2258 const unsigned int n_dofs_per_face = 2259 face_1->get_fe(face_1->nth_active_fe_index(0)).n_dofs_per_face(); 2260 2261 Assert(matrix.m() == 0 || 2262 (first_vector_components.empty() && 2263 matrix.m() == n_dofs_per_face) || 2264 (!first_vector_components.empty() && matrix.m() == spacedim), 2265 ExcMessage( 2266 "The matrix must have either size 0 or spacedim " 2267 "(if first_vector_components is nonempty) " 2268 "or the size must be equal to the # of DoFs on the face " 2269 "(if first_vector_components is empty).")); 2270 } 2271 2272 if (!face_2->has_children()) 2273 { 2274 Assert(face_2->n_active_fe_indices() == 1, ExcInternalError()); 2275 const unsigned int n_dofs_per_face = 2276 face_2->get_fe(face_2->nth_active_fe_index(0)).n_dofs_per_face(); 2277 2278 Assert(matrix.m() == 0 || 2279 (first_vector_components.empty() && 2280 matrix.m() == n_dofs_per_face) || 2281 (!first_vector_components.empty() && matrix.m() == spacedim), 2282 ExcMessage( 2283 "The matrix must have either size 0 or spacedim " 2284 "(if first_vector_components is nonempty) " 2285 "or the size must be equal to the # of DoFs on the face " 2286 "(if first_vector_components is empty).")); 2287 } 2288 #endif 2289 2290 // A lookup table on how to go through the child faces depending on the 2291 // orientation: 2292 2293 static const int lookup_table_2d[2][2] = { 2294 // flip: 2295 {0, 1}, // false 2296 {1, 0}, // true 2297 }; 2298 2299 static const int lookup_table_3d[2][2][2][4] = { 2300 // orientation flip rotation 2301 { 2302 { 2303 {0, 2, 1, 3}, // false false false 2304 {2, 3, 0, 1}, // false false true 2305 }, 2306 { 2307 {3, 1, 2, 0}, // false true false 2308 {1, 0, 3, 2}, // false true true 2309 }, 2310 }, 2311 { 2312 { 2313 {0, 1, 2, 3}, // true false false 2314 {1, 3, 0, 2}, // true false true 2315 }, 2316 { 2317 {3, 2, 1, 0}, // true true false 2318 {2, 0, 3, 1}, // true true true 2319 }, 2320 }, 2321 }; 2322 2323 if (face_1->has_children() && face_2->has_children()) 2324 { 2325 // In the case that both faces have children, we loop over all children 2326 // and apply make_periodicty_constrains recursively: 2327 2328 Assert(face_1->n_children() == 2329 GeometryInfo<dim>::max_children_per_face && 2330 face_2->n_children() == 2331 GeometryInfo<dim>::max_children_per_face, 2332 ExcNotImplemented()); 2333 2334 for (unsigned int i = 0; i < GeometryInfo<dim>::max_children_per_face; 2335 ++i) 2336 { 2337 // Lookup the index for the second face 2338 unsigned int j; 2339 switch (dim) 2340 { 2341 case 2: 2342 j = lookup_table_2d[face_flip][i]; 2343 break; 2344 case 3: 2345 j = lookup_table_3d[face_orientation][face_flip] 2346 [face_rotation][i]; 2347 break; 2348 default: 2349 AssertThrow(false, ExcNotImplemented()); 2350 } 2351 2352 make_periodicity_constraints(face_1->child(i), 2353 face_2->child(j), 2354 affine_constraints, 2355 component_mask, 2356 face_orientation, 2357 face_flip, 2358 face_rotation, 2359 matrix, 2360 first_vector_components, 2361 periodicity_factor); 2362 } 2363 } 2364 else 2365 { 2366 // Otherwise at least one of the two faces is active and we need to do 2367 // some work and enter the constraints! 2368 2369 // The finite element that matters is the one on the active face: 2370 const FiniteElement<dim, spacedim> &fe = 2371 face_1->has_children() ? 2372 face_2->get_fe(face_2->nth_active_fe_index(0)) : 2373 face_1->get_fe(face_1->nth_active_fe_index(0)); 2374 2375 const unsigned int n_dofs_per_face = fe.n_dofs_per_face(); 2376 2377 // Sometimes we just have nothing to do (for all finite elements, or 2378 // systems which accidentally don't have any dofs on the boundary). 2379 if (n_dofs_per_face == 0) 2380 return; 2381 2382 const FullMatrix<double> transformation = 2383 compute_transformation(fe, matrix, first_vector_components); 2384 2385 if (!face_2->has_children()) 2386 { 2387 // Performance hack: We do not need to compute an inverse if the 2388 // matrix is the identity matrix. 2389 if (first_vector_components.empty() && matrix.m() == 0) 2390 { 2391 set_periodicity_constraints(face_2, 2392 face_1, 2393 transformation, 2394 affine_constraints, 2395 component_mask, 2396 face_orientation, 2397 face_flip, 2398 face_rotation, 2399 periodicity_factor); 2400 } 2401 else 2402 { 2403 FullMatrix<double> inverse(transformation.m()); 2404 inverse.invert(transformation); 2405 2406 set_periodicity_constraints(face_2, 2407 face_1, 2408 inverse, 2409 affine_constraints, 2410 component_mask, 2411 face_orientation, 2412 face_flip, 2413 face_rotation, 2414 periodicity_factor); 2415 } 2416 } 2417 else 2418 { 2419 Assert(!face_1->has_children(), ExcInternalError()); 2420 2421 // Important note: 2422 // In 3D we have to take care of the fact that face_rotation gives 2423 // the relative rotation of face_1 to face_2, i.e. we have to invert 2424 // the rotation when constraining face_2 to face_1. Therefore 2425 // face_flip has to be toggled if face_rotation is true: In case of 2426 // inverted orientation, nothing has to be done. 2427 set_periodicity_constraints(face_1, 2428 face_2, 2429 transformation, 2430 affine_constraints, 2431 component_mask, 2432 face_orientation, 2433 face_orientation ? 2434 face_rotation ^ face_flip : 2435 face_flip, 2436 face_rotation, 2437 periodicity_factor); 2438 } 2439 } 2440 } 2441 2442 2443 2444 template <int dim, int spacedim, typename number> 2445 void make_periodicity_constraints(const std::vector<GridTools::PeriodicFacePair<typename DoFHandler<dim,spacedim>::cell_iterator>> & periodic_faces,AffineConstraints<number> & constraints,const ComponentMask & component_mask,const std::vector<unsigned int> & first_vector_components,const number periodicity_factor)2446 make_periodicity_constraints( 2447 const std::vector<GridTools::PeriodicFacePair< 2448 typename DoFHandler<dim, spacedim>::cell_iterator>> &periodic_faces, 2449 AffineConstraints<number> & constraints, 2450 const ComponentMask & component_mask, 2451 const std::vector<unsigned int> &first_vector_components, 2452 const number periodicity_factor) 2453 { 2454 // Loop over all periodic faces... 2455 for (auto &pair : periodic_faces) 2456 { 2457 using FaceIterator = typename DoFHandler<dim, spacedim>::face_iterator; 2458 const FaceIterator face_1 = pair.cell[0]->face(pair.face_idx[0]); 2459 const FaceIterator face_2 = pair.cell[1]->face(pair.face_idx[1]); 2460 2461 Assert(face_1->at_boundary() && face_2->at_boundary(), 2462 ExcInternalError()); 2463 2464 Assert(face_1 != face_2, ExcInternalError()); 2465 2466 // ... and apply the low level make_periodicity_constraints function to 2467 // every matching pair: 2468 make_periodicity_constraints(face_1, 2469 face_2, 2470 constraints, 2471 component_mask, 2472 pair.orientation[0], 2473 pair.orientation[1], 2474 pair.orientation[2], 2475 pair.matrix, 2476 first_vector_components, 2477 periodicity_factor); 2478 } 2479 } 2480 2481 2482 // High level interface variants: 2483 2484 2485 template <int dim, int spacedim, typename number> 2486 void make_periodicity_constraints(const DoFHandler<dim,spacedim> & dof_handler,const types::boundary_id b_id1,const types::boundary_id b_id2,const unsigned int direction,dealii::AffineConstraints<number> & constraints,const ComponentMask & component_mask,const number periodicity_factor)2487 make_periodicity_constraints(const DoFHandler<dim, spacedim> & dof_handler, 2488 const types::boundary_id b_id1, 2489 const types::boundary_id b_id2, 2490 const unsigned int direction, 2491 dealii::AffineConstraints<number> &constraints, 2492 const ComponentMask &component_mask, 2493 const number periodicity_factor) 2494 { 2495 AssertIndexRange(direction, spacedim); 2496 2497 Assert(b_id1 != b_id2, 2498 ExcMessage("The boundary indicators b_id1 and b_id2 must be " 2499 "different to denote different boundaries.")); 2500 2501 std::vector<GridTools::PeriodicFacePair< 2502 typename DoFHandler<dim, spacedim>::cell_iterator>> 2503 matched_faces; 2504 2505 // Collect matching periodic cells on the coarsest level: 2506 GridTools::collect_periodic_faces( 2507 dof_handler, b_id1, b_id2, direction, matched_faces); 2508 2509 make_periodicity_constraints<dim, spacedim>(matched_faces, 2510 constraints, 2511 component_mask, 2512 std::vector<unsigned int>(), 2513 periodicity_factor); 2514 } 2515 2516 2517 2518 template <int dim, int spacedim, typename number> 2519 void make_periodicity_constraints(const DoFHandler<dim,spacedim> & dof_handler,const types::boundary_id b_id,const unsigned int direction,AffineConstraints<number> & constraints,const ComponentMask & component_mask,const number periodicity_factor)2520 make_periodicity_constraints(const DoFHandler<dim, spacedim> &dof_handler, 2521 const types::boundary_id b_id, 2522 const unsigned int direction, 2523 AffineConstraints<number> & constraints, 2524 const ComponentMask & component_mask, 2525 const number periodicity_factor) 2526 { 2527 AssertIndexRange(direction, spacedim); 2528 2529 Assert(dim == spacedim, ExcNotImplemented()); 2530 2531 std::vector<GridTools::PeriodicFacePair< 2532 typename DoFHandler<dim, spacedim>::cell_iterator>> 2533 matched_faces; 2534 2535 // Collect matching periodic cells on the coarsest level: 2536 GridTools::collect_periodic_faces(dof_handler, 2537 b_id, 2538 direction, 2539 matched_faces); 2540 2541 make_periodicity_constraints<dim, spacedim>(matched_faces, 2542 constraints, 2543 component_mask, 2544 std::vector<unsigned int>(), 2545 periodicity_factor); 2546 } 2547 2548 2549 2550 namespace internal 2551 { 2552 namespace Assembler 2553 { 2554 struct Scratch 2555 {}; 2556 2557 2558 template <int dim, int spacedim> 2559 struct CopyData 2560 { 2561 unsigned int dofs_per_cell; 2562 std::vector<types::global_dof_index> parameter_dof_indices; 2563 #ifdef DEAL_II_WITH_MPI 2564 std::vector<dealii::LinearAlgebra::distributed::Vector<double>> 2565 global_parameter_representation; 2566 #else 2567 std::vector<dealii::Vector<double>> global_parameter_representation; 2568 #endif 2569 }; 2570 } // namespace Assembler 2571 2572 namespace 2573 { 2574 /** 2575 * This is a function that is called by the _2 function and that operates 2576 * on one cell only. It is worked in parallel if multhithreading is 2577 * available. 2578 */ 2579 template <int dim, int spacedim> 2580 void compute_intergrid_weights_3(const typename dealii::DoFHandler<dim,spacedim>::active_cell_iterator & cell,const Assembler::Scratch &,Assembler::CopyData<dim,spacedim> & copy_data,const unsigned int coarse_component,const FiniteElement<dim,spacedim> & coarse_fe,const InterGridMap<dealii::DoFHandler<dim,spacedim>> & coarse_to_fine_grid_map,const std::vector<dealii::Vector<double>> & parameter_dofs)2581 compute_intergrid_weights_3( 2582 const typename dealii::DoFHandler<dim, spacedim>::active_cell_iterator 2583 &cell, 2584 const Assembler::Scratch &, 2585 Assembler::CopyData<dim, spacedim> ©_data, 2586 const unsigned int coarse_component, 2587 const FiniteElement<dim, spacedim> &coarse_fe, 2588 const InterGridMap<dealii::DoFHandler<dim, spacedim>> 2589 & coarse_to_fine_grid_map, 2590 const std::vector<dealii::Vector<double>> ¶meter_dofs) 2591 { 2592 // for each cell on the parameter grid: find out which degrees of 2593 // freedom on the fine grid correspond in which way to the degrees of 2594 // freedom on the parameter grid 2595 // 2596 // since for continuous FEs some dofs exist on more than one cell, we 2597 // have to track which ones were already visited. the problem is that if 2598 // we visit a dof first on one cell and compute its weight with respect 2599 // to some global dofs to be non-zero, and later visit the dof again on 2600 // another cell and (since we are on another cell) recompute the weights 2601 // with respect to the same dofs as above to be zero now, we have to 2602 // preserve them. we therefore overwrite all weights if they are nonzero 2603 // and do not enforce zero weights since that might be only due to the 2604 // fact that we are on another cell. 2605 // 2606 // example: 2607 // coarse grid 2608 // | | | 2609 // *-----*-----* 2610 // | cell|cell | 2611 // | 1 | 2 | 2612 // | | | 2613 // 0-----1-----* 2614 // 2615 // fine grid 2616 // | | | | | 2617 // *--*--*--*--* 2618 // | | | | | 2619 // *--*--*--*--* 2620 // | | | | | 2621 // *--x--y--*--* 2622 // 2623 // when on cell 1, we compute the weights of dof 'x' to be 1/2 from 2624 // parameter dofs 0 and 1, respectively. however, when later we are on 2625 // cell 2, we again compute the prolongation of shape function 1 2626 // restricted to cell 2 to the globla grid and find that the weight of 2627 // global dof 'x' now is zero. however, we should not overwrite the old 2628 // value. 2629 // 2630 // we therefore always only set nonzero values. why adding up is not 2631 // useful: dof 'y' would get weight 1 from parameter dof 1 on both cells 2632 // 1 and 2, but the correct weight is nevertheless only 1. 2633 2634 // vector to hold the representation of a single degree of freedom on 2635 // the coarse grid (for the selected fe) on the fine grid 2636 2637 copy_data.dofs_per_cell = coarse_fe.n_dofs_per_cell(); 2638 copy_data.parameter_dof_indices.resize(copy_data.dofs_per_cell); 2639 2640 // get the global indices of the parameter dofs on this parameter grid 2641 // cell 2642 cell->get_dof_indices(copy_data.parameter_dof_indices); 2643 2644 // loop over all dofs on this cell and check whether they are 2645 // interesting for us 2646 for (unsigned int local_dof = 0; local_dof < copy_data.dofs_per_cell; 2647 ++local_dof) 2648 if (coarse_fe.system_to_component_index(local_dof).first == 2649 coarse_component) 2650 { 2651 // the how-many-th parameter is this on this cell? 2652 const unsigned int local_parameter_dof = 2653 coarse_fe.system_to_component_index(local_dof).second; 2654 2655 copy_data.global_parameter_representation[local_parameter_dof] = 2656 0.; 2657 2658 // distribute the representation of @p{local_parameter_dof} on the 2659 // parameter grid cell 2660 // @p{cell} to the global data space 2661 coarse_to_fine_grid_map[cell]->set_dof_values_by_interpolation( 2662 parameter_dofs[local_parameter_dof], 2663 copy_data.global_parameter_representation[local_parameter_dof]); 2664 } 2665 } 2666 2667 2668 2669 /** 2670 * This is a function that is called by the _2 function and that operates 2671 * on one cell only. It is worked in parallel if multhithreading is 2672 * available. 2673 */ 2674 template <int dim, int spacedim> 2675 void copy_intergrid_weights_3(const Assembler::CopyData<dim,spacedim> & copy_data,const unsigned int coarse_component,const FiniteElement<dim,spacedim> & coarse_fe,const std::vector<types::global_dof_index> & weight_mapping,const bool is_called_in_parallel,std::vector<std::map<types::global_dof_index,float>> & weights)2676 copy_intergrid_weights_3( 2677 const Assembler::CopyData<dim, spacedim> & copy_data, 2678 const unsigned int coarse_component, 2679 const FiniteElement<dim, spacedim> & coarse_fe, 2680 const std::vector<types::global_dof_index> &weight_mapping, 2681 const bool is_called_in_parallel, 2682 std::vector<std::map<types::global_dof_index, float>> &weights) 2683 { 2684 unsigned int pos = 0; 2685 for (unsigned int local_dof = 0; local_dof < copy_data.dofs_per_cell; 2686 ++local_dof) 2687 if (coarse_fe.system_to_component_index(local_dof).first == 2688 coarse_component) 2689 { 2690 // now that we've got the global representation of each parameter 2691 // dof, we've only got to clobber the non-zero entries in that 2692 // vector and store the result 2693 // 2694 // what we have learned: if entry @p{i} of the global vector holds 2695 // the value @p{v[i]}, then this is the weight with which the 2696 // present dof contributes to @p{i}. there may be several such 2697 // @p{i}s and their weights' sum should be one. Then, @p{v[i]} 2698 // should be equal to @p{\sum_j w_{ij} p[j]} with @p{p[j]} be the 2699 // values of the degrees of freedom on the coarse grid. we can 2700 // thus compute constraints which link the degrees of freedom 2701 // @p{v[i]} on the fine grid to those on the coarse grid, 2702 // @p{p[j]}. Now to use these as real constraints, rather than as 2703 // additional equations, we have to identify representants among 2704 // the @p{i} for each @p{j}. this will be done by simply taking 2705 // the first @p{i} for which @p{w_{ij}==1}. 2706 // 2707 // guard modification of the weights array by a Mutex. since it 2708 // should happen rather rarely that there are several threads 2709 // operating on different intergrid weights, have only one mutex 2710 // for all of them 2711 for (types::global_dof_index i = 0; 2712 i < copy_data.global_parameter_representation[pos].size(); 2713 ++i) 2714 // set this weight if it belongs to a parameter dof. 2715 if (weight_mapping[i] != numbers::invalid_dof_index) 2716 { 2717 // only overwrite old value if not by zero 2718 if (copy_data.global_parameter_representation[pos](i) != 0) 2719 { 2720 const types::global_dof_index 2721 wi = copy_data.parameter_dof_indices[local_dof], 2722 wj = weight_mapping[i]; 2723 weights[wi][wj] = 2724 copy_data.global_parameter_representation[pos](i); 2725 } 2726 } 2727 else if (!is_called_in_parallel) 2728 { 2729 // Note that when this function operates with distributed 2730 // fine grid, this assertion is switched off since the 2731 // condition does not necessarily hold 2732 Assert(copy_data.global_parameter_representation[pos](i) == 2733 0, 2734 ExcInternalError()); 2735 } 2736 2737 ++pos; 2738 } 2739 } 2740 2741 2742 2743 /** 2744 * This is a helper function that is used in the computation of intergrid 2745 * constraints. See the function for a thorough description of how it 2746 * works. 2747 */ 2748 template <int dim, int spacedim> 2749 void compute_intergrid_weights_2(const dealii::DoFHandler<dim,spacedim> & coarse_grid,const unsigned int coarse_component,const InterGridMap<dealii::DoFHandler<dim,spacedim>> & coarse_to_fine_grid_map,const std::vector<dealii::Vector<double>> & parameter_dofs,const std::vector<types::global_dof_index> & weight_mapping,std::vector<std::map<types::global_dof_index,float>> & weights)2750 compute_intergrid_weights_2( 2751 const dealii::DoFHandler<dim, spacedim> &coarse_grid, 2752 const unsigned int coarse_component, 2753 const InterGridMap<dealii::DoFHandler<dim, spacedim>> 2754 & coarse_to_fine_grid_map, 2755 const std::vector<dealii::Vector<double>> & parameter_dofs, 2756 const std::vector<types::global_dof_index> &weight_mapping, 2757 std::vector<std::map<types::global_dof_index, float>> &weights) 2758 { 2759 Assembler::Scratch scratch; 2760 Assembler::CopyData<dim, spacedim> copy_data; 2761 2762 unsigned int n_interesting_dofs = 0; 2763 for (unsigned int local_dof = 0; 2764 local_dof < coarse_grid.get_fe().n_dofs_per_cell(); 2765 ++local_dof) 2766 if (coarse_grid.get_fe().system_to_component_index(local_dof).first == 2767 coarse_component) 2768 ++n_interesting_dofs; 2769 2770 copy_data.global_parameter_representation.resize(n_interesting_dofs); 2771 2772 bool is_called_in_parallel = false; 2773 for (std::size_t i = 0; 2774 i < copy_data.global_parameter_representation.size(); 2775 ++i) 2776 { 2777 #ifdef DEAL_II_WITH_MPI 2778 MPI_Comm communicator = MPI_COMM_SELF; 2779 try 2780 { 2781 const typename dealii::parallel::TriangulationBase<dim, 2782 spacedim> 2783 &tria = dynamic_cast<const typename dealii::parallel:: 2784 TriangulationBase<dim, spacedim> &>( 2785 coarse_to_fine_grid_map.get_destination_grid() 2786 .get_triangulation()); 2787 communicator = tria.get_communicator(); 2788 is_called_in_parallel = true; 2789 } 2790 catch (std::bad_cast &) 2791 { 2792 // Nothing bad happened: the user used serial Triangulation 2793 } 2794 2795 2796 IndexSet locally_relevant_dofs; 2797 DoFTools::extract_locally_relevant_dofs( 2798 coarse_to_fine_grid_map.get_destination_grid(), 2799 locally_relevant_dofs); 2800 2801 copy_data.global_parameter_representation[i].reinit( 2802 coarse_to_fine_grid_map.get_destination_grid() 2803 .locally_owned_dofs(), 2804 locally_relevant_dofs, 2805 communicator); 2806 #else 2807 const types::global_dof_index n_fine_dofs = weight_mapping.size(); 2808 copy_data.global_parameter_representation[i].reinit(n_fine_dofs); 2809 #endif 2810 } 2811 2812 auto worker = 2813 [coarse_component, 2814 &coarse_grid, 2815 &coarse_to_fine_grid_map, 2816 ¶meter_dofs](const typename dealii::DoFHandler<dim, spacedim>:: 2817 active_cell_iterator & cell, 2818 const Assembler::Scratch & scratch_data, 2819 Assembler::CopyData<dim, spacedim> ©_data) { 2820 compute_intergrid_weights_3<dim, spacedim>(cell, 2821 scratch_data, 2822 copy_data, 2823 coarse_component, 2824 coarse_grid.get_fe(), 2825 coarse_to_fine_grid_map, 2826 parameter_dofs); 2827 }; 2828 2829 auto copier = 2830 [coarse_component, 2831 &coarse_grid, 2832 &weight_mapping, 2833 is_called_in_parallel, 2834 &weights](const Assembler::CopyData<dim, spacedim> ©_data) { 2835 copy_intergrid_weights_3<dim, spacedim>(copy_data, 2836 coarse_component, 2837 coarse_grid.get_fe(), 2838 weight_mapping, 2839 is_called_in_parallel, 2840 weights); 2841 }; 2842 2843 WorkStream::run(coarse_grid.begin_active(), 2844 coarse_grid.end(), 2845 worker, 2846 copier, 2847 scratch, 2848 copy_data); 2849 2850 #ifdef DEAL_II_WITH_MPI 2851 for (std::size_t i = 0; 2852 i < copy_data.global_parameter_representation.size(); 2853 ++i) 2854 copy_data.global_parameter_representation[i].update_ghost_values(); 2855 #endif 2856 } 2857 2858 2859 2860 /** 2861 * This is a helper function that is used in the computation of integrid 2862 * constraints. See the function for a thorough description of how it 2863 * works. 2864 */ 2865 template <int dim, int spacedim> 2866 unsigned int compute_intergrid_weights_1(const dealii::DoFHandler<dim,spacedim> & coarse_grid,const unsigned int coarse_component,const dealii::DoFHandler<dim,spacedim> & fine_grid,const unsigned int fine_component,const InterGridMap<dealii::DoFHandler<dim,spacedim>> & coarse_to_fine_grid_map,std::vector<std::map<types::global_dof_index,float>> & weights,std::vector<types::global_dof_index> & weight_mapping)2867 compute_intergrid_weights_1( 2868 const dealii::DoFHandler<dim, spacedim> &coarse_grid, 2869 const unsigned int coarse_component, 2870 const dealii::DoFHandler<dim, spacedim> &fine_grid, 2871 const unsigned int fine_component, 2872 const InterGridMap<dealii::DoFHandler<dim, spacedim>> 2873 &coarse_to_fine_grid_map, 2874 std::vector<std::map<types::global_dof_index, float>> &weights, 2875 std::vector<types::global_dof_index> & weight_mapping) 2876 { 2877 // aliases to the finite elements used by the dof handlers: 2878 const FiniteElement<dim, spacedim> &coarse_fe = coarse_grid.get_fe(), 2879 &fine_fe = fine_grid.get_fe(); 2880 2881 // global numbers of dofs 2882 const types::global_dof_index n_coarse_dofs = coarse_grid.n_dofs(), 2883 n_fine_dofs = fine_grid.n_dofs(); 2884 2885 // local numbers of dofs 2886 const unsigned int fine_dofs_per_cell = fine_fe.n_dofs_per_cell(); 2887 2888 // alias the number of dofs per cell belonging to the coarse_component 2889 // which is to be the restriction of the fine grid: 2890 const unsigned int coarse_dofs_per_cell_component = 2891 coarse_fe 2892 .base_element( 2893 coarse_fe.component_to_base_index(coarse_component).first) 2894 .n_dofs_per_cell(); 2895 2896 2897 // Try to find out whether the grids stem from the same coarse grid. 2898 // This is a rather crude test, but better than nothing 2899 Assert(coarse_grid.get_triangulation().n_cells(0) == 2900 fine_grid.get_triangulation().n_cells(0), 2901 ExcGridsDontMatch()); 2902 2903 // check whether the map correlates the right objects 2904 Assert(&coarse_to_fine_grid_map.get_source_grid() == &coarse_grid, 2905 ExcGridsDontMatch()); 2906 Assert(&coarse_to_fine_grid_map.get_destination_grid() == &fine_grid, 2907 ExcGridsDontMatch()); 2908 2909 2910 // check whether component numbers are valid 2911 AssertIndexRange(coarse_component, coarse_fe.n_components()); 2912 AssertIndexRange(fine_component, fine_fe.n_components()); 2913 2914 // check whether respective finite elements are equal 2915 Assert(coarse_fe.base_element( 2916 coarse_fe.component_to_base_index(coarse_component).first) == 2917 fine_fe.base_element( 2918 fine_fe.component_to_base_index(fine_component).first), 2919 ExcFiniteElementsDontMatch()); 2920 2921 #ifdef DEBUG 2922 // if in debug mode, check whether the coarse grid is indeed coarser 2923 // everywhere than the fine grid 2924 for (const auto &cell : coarse_grid.active_cell_iterators()) 2925 Assert(cell->level() <= coarse_to_fine_grid_map[cell]->level(), 2926 ExcGridNotCoarser()); 2927 #endif 2928 2929 /* 2930 * From here on: the term `parameter' refers to the selected component 2931 * on the coarse grid and its analogon on the fine grid. The naming of 2932 * variables containing this term is due to the fact that 2933 * `selected_component' is longer, but also due to the fact that the 2934 * code of this function was initially written for a program where the 2935 * component which we wanted to match between grids was actually the 2936 * `parameter' variable. 2937 * 2938 * Likewise, the terms `parameter grid' and `state grid' refer to the 2939 * coarse and fine grids, respectively. 2940 * 2941 * Changing the names of variables would in principle be a good idea, 2942 * but would not make things simpler and would be another source of 2943 * errors. If anyone feels like doing so: patches would be welcome! 2944 */ 2945 2946 2947 2948 // set up vectors of cell-local data; each vector represents one degree 2949 // of freedom of the coarse-grid variable in the fine-grid element 2950 std::vector<dealii::Vector<double>> parameter_dofs( 2951 coarse_dofs_per_cell_component, 2952 dealii::Vector<double>(fine_dofs_per_cell)); 2953 // for each coarse dof: find its position within the fine element and 2954 // set this value to one in the respective vector (all other values are 2955 // zero by construction) 2956 for (unsigned int local_coarse_dof = 0; 2957 local_coarse_dof < coarse_dofs_per_cell_component; 2958 ++local_coarse_dof) 2959 for (unsigned int fine_dof = 0; fine_dof < fine_fe.n_dofs_per_cell(); 2960 ++fine_dof) 2961 if (fine_fe.system_to_component_index(fine_dof) == 2962 std::make_pair(fine_component, local_coarse_dof)) 2963 { 2964 parameter_dofs[local_coarse_dof](fine_dof) = 1.; 2965 break; 2966 } 2967 2968 2969 // find out how many DoFs there are on the grids belonging to the 2970 // components we want to match 2971 unsigned int n_parameters_on_fine_grid = 0; 2972 { 2973 // have a flag for each dof on the fine grid and set it to true if 2974 // this is an interesting dof. finally count how many true's there 2975 std::vector<bool> dof_is_interesting(fine_grid.n_dofs(), false); 2976 std::vector<types::global_dof_index> local_dof_indices( 2977 fine_fe.n_dofs_per_cell()); 2978 2979 for (const auto &cell : fine_grid.active_cell_iterators()) 2980 if (cell->is_locally_owned()) 2981 { 2982 cell->get_dof_indices(local_dof_indices); 2983 for (unsigned int i = 0; i < fine_fe.n_dofs_per_cell(); ++i) 2984 if (fine_fe.system_to_component_index(i).first == 2985 fine_component) 2986 dof_is_interesting[local_dof_indices[i]] = true; 2987 } 2988 2989 n_parameters_on_fine_grid = std::count(dof_is_interesting.begin(), 2990 dof_is_interesting.end(), 2991 true); 2992 } 2993 2994 2995 // set up the weights mapping 2996 weights.clear(); 2997 weights.resize(n_coarse_dofs); 2998 2999 weight_mapping.clear(); 3000 weight_mapping.resize(n_fine_dofs, numbers::invalid_dof_index); 3001 3002 { 3003 std::vector<types::global_dof_index> local_dof_indices( 3004 fine_fe.n_dofs_per_cell()); 3005 unsigned int next_free_index = 0; 3006 for (const auto &cell : fine_grid.active_cell_iterators()) 3007 if (cell->is_locally_owned()) 3008 { 3009 cell->get_dof_indices(local_dof_indices); 3010 for (unsigned int i = 0; i < fine_fe.n_dofs_per_cell(); ++i) 3011 // if this DoF is a parameter dof and has not yet been 3012 // numbered, then do so 3013 if ((fine_fe.system_to_component_index(i).first == 3014 fine_component) && 3015 (weight_mapping[local_dof_indices[i]] == 3016 numbers::invalid_dof_index)) 3017 { 3018 weight_mapping[local_dof_indices[i]] = next_free_index; 3019 ++next_free_index; 3020 } 3021 } 3022 3023 Assert(next_free_index == n_parameters_on_fine_grid, 3024 ExcInternalError()); 3025 } 3026 3027 3028 // for each cell on the parameter grid: find out which degrees of 3029 // freedom on the fine grid correspond in which way to the degrees of 3030 // freedom on the parameter grid 3031 // 3032 // do this in a separate function to allow for multithreading there. see 3033 // this function also if you want to read more information on the 3034 // algorithm used. 3035 compute_intergrid_weights_2(coarse_grid, 3036 coarse_component, 3037 coarse_to_fine_grid_map, 3038 parameter_dofs, 3039 weight_mapping, 3040 weights); 3041 3042 3043 // ok, now we have all weights for each dof on the fine grid. if in 3044 // debug mode lets see if everything went smooth, i.e. each dof has sum 3045 // of weights one 3046 // 3047 // in other words this means that if the sum of all shape functions on 3048 // the parameter grid is one (which is always the case), then the 3049 // representation on the state grid should be as well (division of 3050 // unity) 3051 // 3052 // if the parameter grid has more than one component, then the 3053 // respective dofs of the other components have sum of weights zero, of 3054 // course. we do not explicitly ask which component a dof belongs to, 3055 // but this at least tests some errors 3056 #ifdef DEBUG 3057 for (unsigned int col = 0; col < n_parameters_on_fine_grid; ++col) 3058 { 3059 double sum = 0; 3060 for (types::global_dof_index row = 0; row < n_coarse_dofs; ++row) 3061 if (weights[row].find(col) != weights[row].end()) 3062 sum += weights[row][col]; 3063 Assert((std::fabs(sum - 1) < 1.e-12) || 3064 ((coarse_fe.n_components() > 1) && (sum == 0)), 3065 ExcInternalError()); 3066 } 3067 #endif 3068 3069 3070 return n_parameters_on_fine_grid; 3071 } 3072 3073 3074 } // namespace 3075 } // namespace internal 3076 3077 3078 3079 template <int dim, int spacedim> 3080 void compute_intergrid_constraints(const DoFHandler<dim,spacedim> & coarse_grid,const unsigned int coarse_component,const DoFHandler<dim,spacedim> & fine_grid,const unsigned int fine_component,const InterGridMap<DoFHandler<dim,spacedim>> & coarse_to_fine_grid_map,AffineConstraints<double> & constraints)3081 compute_intergrid_constraints( 3082 const DoFHandler<dim, spacedim> & coarse_grid, 3083 const unsigned int coarse_component, 3084 const DoFHandler<dim, spacedim> & fine_grid, 3085 const unsigned int fine_component, 3086 const InterGridMap<DoFHandler<dim, spacedim>> &coarse_to_fine_grid_map, 3087 AffineConstraints<double> & constraints) 3088 { 3089 Assert(coarse_grid.get_fe_collection().size() == 1 && 3090 fine_grid.get_fe_collection().size() == 1, 3091 ExcMessage("This function is not yet implemented for DoFHandlers " 3092 "using hp capabilities.")); 3093 // store the weights with which a dof on the parameter grid contributes to a 3094 // dof on the fine grid. see the long doc below for more info 3095 // 3096 // allocate as many rows as there are parameter dofs on the coarse grid and 3097 // as many columns as there are parameter dofs on the fine grid. 3098 // 3099 // weight_mapping is used to map the global (fine grid) parameter dof 3100 // indices to the columns 3101 // 3102 // in the original implementation, the weights array was actually of 3103 // FullMatrix<double> type. this wasted huge amounts of memory, but was 3104 // fast. nonetheless, since the memory consumption was quadratic in the 3105 // number of degrees of freedom, this was not very practical, so we now use 3106 // a vector of rows of the matrix, and in each row a vector of pairs 3107 // (colnum,value). this seems like the best tradeoff between memory and 3108 // speed, as it is now linear in memory and still fast enough. 3109 // 3110 // to save some memory and since the weights are usually (negative) powers 3111 // of 2, we choose the value type of the matrix to be @p{float} rather than 3112 // @p{double}. 3113 std::vector<std::map<types::global_dof_index, float>> weights; 3114 3115 // this is this mapping. there is one entry for each dof on the fine grid; 3116 // if it is a parameter dof, then its value is the column in weights for 3117 // that parameter dof, if it is any other dof, then its value is -1, 3118 // indicating an error 3119 std::vector<types::global_dof_index> weight_mapping; 3120 3121 const unsigned int n_parameters_on_fine_grid = 3122 internal::compute_intergrid_weights_1(coarse_grid, 3123 coarse_component, 3124 fine_grid, 3125 fine_component, 3126 coarse_to_fine_grid_map, 3127 weights, 3128 weight_mapping); 3129 (void)n_parameters_on_fine_grid; 3130 3131 // global numbers of dofs 3132 const types::global_dof_index n_coarse_dofs = coarse_grid.n_dofs(), 3133 n_fine_dofs = fine_grid.n_dofs(); 3134 3135 3136 // get an array in which we store which dof on the coarse grid is a 3137 // parameter and which is not 3138 IndexSet coarse_dof_is_parameter; 3139 { 3140 std::vector<bool> mask(coarse_grid.get_fe(0).n_components(), false); 3141 mask[coarse_component] = true; 3142 3143 coarse_dof_is_parameter = 3144 extract_dofs<dim, spacedim>(coarse_grid, ComponentMask(mask)); 3145 } 3146 3147 // now we know that the weights in each row constitute a constraint. enter 3148 // this into the constraints object 3149 // 3150 // first task: for each parameter dof on the parameter grid, find a 3151 // representant on the fine, global grid. this is possible since we use 3152 // conforming finite element. we take this representant to be the first 3153 // element in this row with weight identical to one. the representant will 3154 // become an unconstrained degree of freedom, while all others will be 3155 // constrained to this dof (and possibly others) 3156 std::vector<types::global_dof_index> representants( 3157 n_coarse_dofs, numbers::invalid_dof_index); 3158 for (types::global_dof_index parameter_dof = 0; 3159 parameter_dof < n_coarse_dofs; 3160 ++parameter_dof) 3161 if (coarse_dof_is_parameter.is_element(parameter_dof)) 3162 { 3163 // if this is the line of a parameter dof on the coarse grid, then it 3164 // should have at least one dependent node on the fine grid 3165 Assert(weights[parameter_dof].size() > 0, ExcInternalError()); 3166 3167 // find the column where the representant is mentioned 3168 std::map<types::global_dof_index, float>::const_iterator i = 3169 weights[parameter_dof].begin(); 3170 for (; i != weights[parameter_dof].end(); ++i) 3171 if (i->second == 1) 3172 break; 3173 Assert(i != weights[parameter_dof].end(), ExcInternalError()); 3174 const types::global_dof_index column = i->first; 3175 3176 // now we know in which column of weights the representant is, but we 3177 // don't know its global index. get it using the inverse operation of 3178 // the weight_mapping 3179 types::global_dof_index global_dof = 0; 3180 for (; global_dof < weight_mapping.size(); ++global_dof) 3181 if (weight_mapping[global_dof] == 3182 static_cast<types::global_dof_index>(column)) 3183 break; 3184 Assert(global_dof < weight_mapping.size(), ExcInternalError()); 3185 3186 // now enter the representants global index into our list 3187 representants[parameter_dof] = global_dof; 3188 } 3189 else 3190 { 3191 // consistency check: if this is no parameter dof on the coarse grid, 3192 // then the respective row must be empty! 3193 Assert(weights[parameter_dof].size() == 0, ExcInternalError()); 3194 } 3195 3196 3197 3198 // note for people that want to optimize this function: the largest part of 3199 // the computing time is spent in the following, rather innocent block of 3200 // code. basically, it must be the AffineConstraints::add_entry call which 3201 // takes the bulk of the time, but it is not known to the author how to make 3202 // it faster... 3203 std::vector<std::pair<types::global_dof_index, double>> constraint_line; 3204 for (types::global_dof_index global_dof = 0; global_dof < n_fine_dofs; 3205 ++global_dof) 3206 if (weight_mapping[global_dof] != numbers::invalid_dof_index) 3207 // this global dof is a parameter dof, so it may carry a constraint note 3208 // that for each global dof, the sum of weights shall be one, so we can 3209 // find out whether this dof is constrained in the following way: if the 3210 // only weight in this row is a one, and the representant for the 3211 // parameter dof of the line in which this one is is the present dof, 3212 // then we consider this dof to be unconstrained. otherwise, all other 3213 // dofs are constrained 3214 { 3215 const types::global_dof_index col = weight_mapping[global_dof]; 3216 Assert(col < n_parameters_on_fine_grid, ExcInternalError()); 3217 3218 types::global_dof_index first_used_row = 0; 3219 3220 { 3221 Assert(weights.size() > 0, ExcInternalError()); 3222 std::map<types::global_dof_index, float>::const_iterator col_entry = 3223 weights[0].end(); 3224 for (; first_used_row < n_coarse_dofs; ++first_used_row) 3225 { 3226 col_entry = weights[first_used_row].find(col); 3227 if (col_entry != weights[first_used_row].end()) 3228 break; 3229 } 3230 3231 Assert(col_entry != weights[first_used_row].end(), 3232 ExcInternalError()); 3233 3234 if ((col_entry->second == 1) && 3235 (representants[first_used_row] == global_dof)) 3236 // dof unconstrained or constrained to itself (in case this cell 3237 // is mapped to itself, rather than to children of itself) 3238 continue; 3239 } 3240 3241 3242 // otherwise enter all constraints 3243 constraints.add_line(global_dof); 3244 3245 constraint_line.clear(); 3246 for (types::global_dof_index row = first_used_row; 3247 row < n_coarse_dofs; 3248 ++row) 3249 { 3250 const std::map<types::global_dof_index, float>::const_iterator j = 3251 weights[row].find(col); 3252 if ((j != weights[row].end()) && (j->second != 0)) 3253 constraint_line.emplace_back(representants[row], j->second); 3254 } 3255 3256 constraints.add_entries(global_dof, constraint_line); 3257 } 3258 } 3259 3260 3261 3262 template <int dim, int spacedim> 3263 void compute_intergrid_transfer_representation(const DoFHandler<dim,spacedim> & coarse_grid,const unsigned int coarse_component,const DoFHandler<dim,spacedim> & fine_grid,const unsigned int fine_component,const InterGridMap<DoFHandler<dim,spacedim>> & coarse_to_fine_grid_map,std::vector<std::map<types::global_dof_index,float>> & transfer_representation)3264 compute_intergrid_transfer_representation( 3265 const DoFHandler<dim, spacedim> & coarse_grid, 3266 const unsigned int coarse_component, 3267 const DoFHandler<dim, spacedim> & fine_grid, 3268 const unsigned int fine_component, 3269 const InterGridMap<DoFHandler<dim, spacedim>> &coarse_to_fine_grid_map, 3270 std::vector<std::map<types::global_dof_index, float>> 3271 &transfer_representation) 3272 { 3273 Assert(coarse_grid.get_fe_collection().size() == 1 && 3274 fine_grid.get_fe_collection().size() == 1, 3275 ExcMessage("This function is not yet implemented for DoFHandlers " 3276 "using hp capabilities.")); 3277 // store the weights with which a dof on the parameter grid contributes to a 3278 // dof on the fine grid. see the long doc below for more info 3279 // 3280 // allocate as many rows as there are parameter dofs on the coarse grid and 3281 // as many columns as there are parameter dofs on the fine grid. 3282 // 3283 // weight_mapping is used to map the global (fine grid) parameter dof 3284 // indices to the columns 3285 // 3286 // in the original implementation, the weights array was actually of 3287 // FullMatrix<double> type. this wasted huge amounts of memory, but was 3288 // fast. nonetheless, since the memory consumption was quadratic in the 3289 // number of degrees of freedom, this was not very practical, so we now use 3290 // a vector of rows of the matrix, and in each row a vector of pairs 3291 // (colnum,value). this seems like the best tradeoff between memory and 3292 // speed, as it is now linear in memory and still fast enough. 3293 // 3294 // to save some memory and since the weights are usually (negative) powers 3295 // of 2, we choose the value type of the matrix to be @p{float} rather than 3296 // @p{double}. 3297 std::vector<std::map<types::global_dof_index, float>> weights; 3298 3299 // this is this mapping. there is one entry for each dof on the fine grid; 3300 // if it is a parameter dof, then its value is the column in weights for 3301 // that parameter dof, if it is any other dof, then its value is -1, 3302 // indicating an error 3303 std::vector<types::global_dof_index> weight_mapping; 3304 3305 internal::compute_intergrid_weights_1(coarse_grid, 3306 coarse_component, 3307 fine_grid, 3308 fine_component, 3309 coarse_to_fine_grid_map, 3310 weights, 3311 weight_mapping); 3312 3313 // now compute the requested representation 3314 const types::global_dof_index n_global_parm_dofs = 3315 std::count_if(weight_mapping.begin(), 3316 weight_mapping.end(), 3317 [](const types::global_dof_index dof) { 3318 return dof != numbers::invalid_dof_index; 3319 }); 3320 3321 // first construct the inverse mapping of weight_mapping 3322 std::vector<types::global_dof_index> inverse_weight_mapping( 3323 n_global_parm_dofs, numbers::invalid_dof_index); 3324 for (types::global_dof_index i = 0; i < weight_mapping.size(); ++i) 3325 { 3326 const types::global_dof_index parameter_dof = weight_mapping[i]; 3327 // if this global dof is a parameter 3328 if (parameter_dof != numbers::invalid_dof_index) 3329 { 3330 Assert(parameter_dof < n_global_parm_dofs, ExcInternalError()); 3331 Assert((inverse_weight_mapping[parameter_dof] == 3332 numbers::invalid_dof_index), 3333 ExcInternalError()); 3334 3335 inverse_weight_mapping[parameter_dof] = i; 3336 } 3337 } 3338 3339 // next copy over weights array and replace respective numbers 3340 const types::global_dof_index n_rows = weight_mapping.size(); 3341 3342 transfer_representation.clear(); 3343 transfer_representation.resize(n_rows); 3344 3345 const types::global_dof_index n_coarse_dofs = coarse_grid.n_dofs(); 3346 for (types::global_dof_index i = 0; i < n_coarse_dofs; ++i) 3347 { 3348 std::map<types::global_dof_index, float>::const_iterator j = 3349 weights[i].begin(); 3350 for (; j != weights[i].end(); ++j) 3351 { 3352 const types::global_dof_index p = inverse_weight_mapping[j->first]; 3353 Assert(p < n_rows, ExcInternalError()); 3354 3355 transfer_representation[p][i] = j->second; 3356 } 3357 } 3358 } 3359 3360 3361 3362 template <int dim, int spacedim, typename number> 3363 void make_zero_boundary_constraints(const DoFHandler<dim,spacedim> & dof,const types::boundary_id boundary_id,AffineConstraints<number> & zero_boundary_constraints,const ComponentMask & component_mask)3364 make_zero_boundary_constraints( 3365 const DoFHandler<dim, spacedim> &dof, 3366 const types::boundary_id boundary_id, 3367 AffineConstraints<number> & zero_boundary_constraints, 3368 const ComponentMask & component_mask) 3369 { 3370 Assert(component_mask.represents_n_components(dof.get_fe(0).n_components()), 3371 ExcMessage("The number of components in the mask has to be either " 3372 "zero or equal to the number of components in the finite " 3373 "element.")); 3374 3375 const unsigned int n_components = dof.get_fe_collection().n_components(); 3376 3377 Assert(component_mask.n_selected_components(n_components) > 0, 3378 ComponentMask::ExcNoComponentSelected()); 3379 3380 // a field to store the indices on the face 3381 std::vector<types::global_dof_index> face_dofs; 3382 face_dofs.reserve(dof.get_fe_collection().max_dofs_per_face()); 3383 // a field to store the indices on the cell 3384 std::vector<types::global_dof_index> cell_dofs; 3385 cell_dofs.reserve(dof.get_fe_collection().max_dofs_per_cell()); 3386 3387 typename DoFHandler<dim, spacedim>::active_cell_iterator 3388 cell = dof.begin_active(), 3389 endc = dof.end(); 3390 for (; cell != endc; ++cell) 3391 if (!cell->is_artificial() && cell->at_boundary()) 3392 { 3393 const FiniteElement<dim, spacedim> &fe = cell->get_fe(); 3394 3395 // get global indices of dofs on the cell 3396 cell_dofs.resize(fe.n_dofs_per_cell()); 3397 cell->get_dof_indices(cell_dofs); 3398 3399 for (const auto face_no : cell->face_indices()) 3400 { 3401 const typename DoFHandler<dim, spacedim>::face_iterator face = 3402 cell->face(face_no); 3403 3404 // if face is on the boundary and satisfies the correct boundary 3405 // id property 3406 if (face->at_boundary() && 3407 ((boundary_id == numbers::invalid_boundary_id) || 3408 (face->boundary_id() == boundary_id))) 3409 { 3410 // get indices and physical location on this face 3411 face_dofs.resize(fe.n_dofs_per_face()); 3412 face->get_dof_indices(face_dofs, cell->active_fe_index()); 3413 3414 // enter those dofs into the list that match the component 3415 // signature. 3416 for (const types::global_dof_index face_dof : face_dofs) 3417 { 3418 // Find out if a dof has a contribution in this component, 3419 // and if so, add it to the list 3420 const std::vector<types::global_dof_index>::iterator 3421 it_index_on_cell = std::find(cell_dofs.begin(), 3422 cell_dofs.end(), 3423 face_dof); 3424 Assert(it_index_on_cell != cell_dofs.end(), 3425 ExcInvalidIterator()); 3426 const unsigned int index_on_cell = 3427 std::distance(cell_dofs.begin(), it_index_on_cell); 3428 const ComponentMask &nonzero_component_array = 3429 cell->get_fe().get_nonzero_components(index_on_cell); 3430 bool nonzero = false; 3431 for (unsigned int c = 0; c < n_components; ++c) 3432 if (nonzero_component_array[c] && component_mask[c]) 3433 { 3434 nonzero = true; 3435 break; 3436 } 3437 3438 if (nonzero) 3439 zero_boundary_constraints.add_line(face_dof); 3440 } 3441 } 3442 } 3443 } 3444 } 3445 3446 3447 3448 template <int dim, int spacedim, typename number> 3449 void make_zero_boundary_constraints(const DoFHandler<dim,spacedim> & dof,AffineConstraints<number> & zero_boundary_constraints,const ComponentMask & component_mask)3450 make_zero_boundary_constraints( 3451 const DoFHandler<dim, spacedim> &dof, 3452 AffineConstraints<number> & zero_boundary_constraints, 3453 const ComponentMask & component_mask) 3454 { 3455 make_zero_boundary_constraints(dof, 3456 numbers::invalid_boundary_id, 3457 zero_boundary_constraints, 3458 component_mask); 3459 } 3460 3461 3462 } // end of namespace DoFTools 3463 3464 3465 3466 // explicit instantiations 3467 3468 #include "dof_tools_constraints.inst" 3469 3470 3471 3472 DEAL_II_NAMESPACE_CLOSE 3473