1 // --------------------------------------------------------------------- 2 // 3 // Copyright (C) 2000 - 2020 by the deal.II authors 4 // 5 // This file is part of the deal.II library. 6 // 7 // The deal.II library is free software; you can use it, redistribute 8 // it, and/or modify it under the terms of the GNU Lesser General 9 // Public License as published by the Free Software Foundation; either 10 // version 2.1 of the License, or (at your option) any later version. 11 // The full text of the license can be found in the file LICENSE.md at 12 // the top level directory of deal.II. 13 // 14 // --------------------------------------------------------------------- 15 16 #ifndef dealii_fe_tools_templates_H 17 #define dealii_fe_tools_templates_H 18 19 20 #include <deal.II/base/config.h> 21 22 #include <deal.II/base/index_set.h> 23 #include <deal.II/base/qprojector.h> 24 #include <deal.II/base/quadrature_lib.h> 25 #include <deal.II/base/thread_management.h> 26 #include <deal.II/base/utilities.h> 27 28 #include <deal.II/dofs/dof_accessor.h> 29 #include <deal.II/dofs/dof_handler.h> 30 #include <deal.II/dofs/dof_tools.h> 31 32 #include <deal.II/fe/fe.h> 33 #include <deal.II/fe/fe_abf.h> 34 #include <deal.II/fe/fe_bdm.h> 35 #include <deal.II/fe/fe_bernstein.h> 36 #include <deal.II/fe/fe_dg_vector.h> 37 #include <deal.II/fe/fe_dgp.h> 38 #include <deal.II/fe/fe_dgp_monomial.h> 39 #include <deal.II/fe/fe_dgp_nonparametric.h> 40 #include <deal.II/fe/fe_dgq.h> 41 #include <deal.II/fe/fe_face.h> 42 #include <deal.II/fe/fe_nedelec.h> 43 #include <deal.II/fe/fe_nedelec_sz.h> 44 #include <deal.II/fe/fe_nothing.h> 45 #include <deal.II/fe/fe_q.h> 46 #include <deal.II/fe/fe_q_bubbles.h> 47 #include <deal.II/fe/fe_q_dg0.h> 48 #include <deal.II/fe/fe_q_hierarchical.h> 49 #include <deal.II/fe/fe_q_iso_q1.h> 50 #include <deal.II/fe/fe_rannacher_turek.h> 51 #include <deal.II/fe/fe_raviart_thomas.h> 52 #include <deal.II/fe/fe_rt_bubbles.h> 53 #include <deal.II/fe/fe_system.h> 54 #include <deal.II/fe/fe_tools.h> 55 #include <deal.II/fe/fe_values.h> 56 #include <deal.II/fe/mapping_cartesian.h> 57 #include <deal.II/fe/mapping_q1.h> 58 59 #include <deal.II/grid/grid_generator.h> 60 #include <deal.II/grid/tria.h> 61 #include <deal.II/grid/tria_iterator.h> 62 63 #include <deal.II/hp/dof_handler.h> 64 65 #include <deal.II/lac/full_matrix.h> 66 #include <deal.II/lac/householder.h> 67 68 #include <cctype> 69 #include <iostream> 70 #include <memory> 71 72 73 DEAL_II_NAMESPACE_OPEN 74 75 namespace FETools 76 { 77 namespace Compositing 78 { 79 template <int dim, int spacedim> 80 FiniteElementData<dim> multiply_dof_numbers(const std::vector<const FiniteElement<dim,spacedim> * > & fes,const std::vector<unsigned int> & multiplicities,const bool do_tensor_product)81 multiply_dof_numbers( 82 const std::vector<const FiniteElement<dim, spacedim> *> &fes, 83 const std::vector<unsigned int> & multiplicities, 84 const bool do_tensor_product) 85 { 86 AssertDimension(fes.size(), multiplicities.size()); 87 88 unsigned int multiplied_dofs_per_vertex = 0; 89 unsigned int multiplied_dofs_per_line = 0; 90 unsigned int multiplied_dofs_per_quad = 0; 91 unsigned int multiplied_dofs_per_hex = 0; 92 93 unsigned int multiplied_n_components = 0; 94 95 unsigned int degree = 0; // degree is the maximal degree of the components 96 97 unsigned int n_components = 0; 98 // Get the number of components from the first given finite element. 99 for (unsigned int i = 0; i < fes.size(); i++) 100 if (multiplicities[i] > 0) 101 { 102 n_components = fes[i]->n_components(); 103 break; 104 } 105 106 for (unsigned int i = 0; i < fes.size(); i++) 107 if (multiplicities[i] > 0) 108 { 109 // TODO: the implementation makes the assumption that all faces have 110 // the same number of dofs -> don't construct DPO but 111 // PrecomputedFiniteElementData 112 AssertDimension(fes[i]->n_unique_quads(), 1); 113 114 multiplied_dofs_per_vertex += 115 fes[i]->n_dofs_per_vertex() * multiplicities[i]; 116 multiplied_dofs_per_line += 117 fes[i]->n_dofs_per_line() * multiplicities[i]; 118 multiplied_dofs_per_quad += 119 fes[i]->n_dofs_per_quad(0) * multiplicities[i]; 120 multiplied_dofs_per_hex += 121 fes[i]->n_dofs_per_hex() * multiplicities[i]; 122 123 multiplied_n_components += 124 fes[i]->n_components() * multiplicities[i]; 125 126 Assert(do_tensor_product || 127 (n_components == fes[i]->n_components()), 128 ExcDimensionMismatch(n_components, fes[i]->n_components())); 129 130 degree = std::max(degree, fes[i]->tensor_degree()); 131 } 132 133 // assume conformity of the first finite element and then take away 134 // bits as indicated by the base elements. if all multiplicities 135 // happen to be zero, then it doesn't matter what we set it to. 136 typename FiniteElementData<dim>::Conformity total_conformity = 137 typename FiniteElementData<dim>::Conformity(); 138 { 139 unsigned int index = 0; 140 for (index = 0; index < fes.size(); ++index) 141 if (multiplicities[index] > 0) 142 { 143 total_conformity = fes[index]->conforming_space; 144 break; 145 } 146 147 for (; index < fes.size(); ++index) 148 if (multiplicities[index] > 0) 149 total_conformity = typename FiniteElementData<dim>::Conformity( 150 total_conformity & fes[index]->conforming_space); 151 } 152 153 std::vector<unsigned int> dpo; 154 dpo.push_back(multiplied_dofs_per_vertex); 155 dpo.push_back(multiplied_dofs_per_line); 156 if (dim > 1) 157 dpo.push_back(multiplied_dofs_per_quad); 158 if (dim > 2) 159 dpo.push_back(multiplied_dofs_per_hex); 160 161 BlockIndices block_indices(0, 0); 162 163 for (unsigned int base = 0; base < fes.size(); ++base) 164 for (unsigned int m = 0; m < multiplicities[base]; ++m) 165 block_indices.push_back(fes[base]->n_dofs_per_cell()); 166 167 return FiniteElementData<dim>(dpo, 168 fes.front()->reference_cell_type(), 169 (do_tensor_product ? 170 multiplied_n_components : 171 n_components), 172 degree, 173 total_conformity, 174 block_indices); 175 } 176 177 178 179 template <int dim, int spacedim> 180 FiniteElementData<dim> multiply_dof_numbers(const FiniteElement<dim,spacedim> * fe1,const unsigned int N1,const FiniteElement<dim,spacedim> * fe2,const unsigned int N2,const FiniteElement<dim,spacedim> * fe3,const unsigned int N3,const FiniteElement<dim,spacedim> * fe4,const unsigned int N4,const FiniteElement<dim,spacedim> * fe5,const unsigned int N5)181 multiply_dof_numbers(const FiniteElement<dim, spacedim> *fe1, 182 const unsigned int N1, 183 const FiniteElement<dim, spacedim> *fe2, 184 const unsigned int N2, 185 const FiniteElement<dim, spacedim> *fe3, 186 const unsigned int N3, 187 const FiniteElement<dim, spacedim> *fe4, 188 const unsigned int N4, 189 const FiniteElement<dim, spacedim> *fe5, 190 const unsigned int N5) 191 { 192 std::vector<const FiniteElement<dim, spacedim> *> fes; 193 fes.push_back(fe1); 194 fes.push_back(fe2); 195 fes.push_back(fe3); 196 fes.push_back(fe4); 197 fes.push_back(fe5); 198 199 std::vector<unsigned int> mult; 200 mult.push_back(N1); 201 mult.push_back(N2); 202 mult.push_back(N3); 203 mult.push_back(N4); 204 mult.push_back(N5); 205 return multiply_dof_numbers(fes, mult); 206 } 207 208 209 210 template <int dim, int spacedim> 211 std::vector<bool> compute_restriction_is_additive_flags(const std::vector<const FiniteElement<dim,spacedim> * > & fes,const std::vector<unsigned int> & multiplicities)212 compute_restriction_is_additive_flags( 213 const std::vector<const FiniteElement<dim, spacedim> *> &fes, 214 const std::vector<unsigned int> & multiplicities) 215 { 216 AssertDimension(fes.size(), multiplicities.size()); 217 218 // first count the number of dofs and components that will emerge from the 219 // given FEs 220 unsigned int n_shape_functions = 0; 221 for (unsigned int i = 0; i < fes.size(); ++i) 222 if (multiplicities[i] > 0) // check needed as fe might be nullptr 223 n_shape_functions += fes[i]->n_dofs_per_cell() * multiplicities[i]; 224 225 // generate the array that will hold the output 226 std::vector<bool> retval(n_shape_functions, false); 227 228 // finally go through all the shape functions of the base elements, and 229 // copy their flags. this somehow copies the code in build_cell_table, 230 // which is not nice as it uses too much implicit knowledge about the 231 // layout of the individual bases in the composed FE, but there seems no 232 // way around... 233 // 234 // for each shape function, copy the flags from the base element to this 235 // one, taking into account multiplicities, and other complications 236 unsigned int total_index = 0; 237 for (const unsigned int vertex_number : 238 ReferenceCell::internal::Info::get_cell( 239 fes.front()->reference_cell_type()) 240 .vertex_indices()) 241 { 242 for (unsigned int base = 0; base < fes.size(); ++base) 243 for (unsigned int m = 0; m < multiplicities[base]; ++m) 244 for (unsigned int local_index = 0; 245 local_index < fes[base]->n_dofs_per_vertex(); 246 ++local_index, ++total_index) 247 { 248 const unsigned int index_in_base = 249 (fes[base]->n_dofs_per_vertex() * vertex_number + 250 local_index); 251 252 Assert(index_in_base < fes[base]->n_dofs_per_cell(), 253 ExcInternalError()); 254 retval[total_index] = 255 fes[base]->restriction_is_additive(index_in_base); 256 } 257 } 258 259 // 2. Lines 260 for (const unsigned int line_number : 261 ReferenceCell::internal::Info::get_cell( 262 fes.front()->reference_cell_type()) 263 .line_indices()) 264 { 265 for (unsigned int base = 0; base < fes.size(); ++base) 266 for (unsigned int m = 0; m < multiplicities[base]; ++m) 267 for (unsigned int local_index = 0; 268 local_index < fes[base]->n_dofs_per_line(); 269 ++local_index, ++total_index) 270 { 271 const unsigned int index_in_base = 272 (fes[base]->n_dofs_per_line() * line_number + local_index + 273 fes[base]->get_first_line_index()); 274 275 Assert(index_in_base < fes[base]->n_dofs_per_cell(), 276 ExcInternalError()); 277 retval[total_index] = 278 fes[base]->restriction_is_additive(index_in_base); 279 } 280 } 281 282 // 3. Quads 283 for (unsigned int quad_number = 0; 284 quad_number < (dim == 2 ? 285 1 : 286 (dim == 3 ? ReferenceCell::internal::Info::get_cell( 287 fes.front()->reference_cell_type()) 288 .n_faces() : 289 0)); 290 ++quad_number) 291 { 292 for (unsigned int base = 0; base < fes.size(); ++base) 293 for (unsigned int m = 0; m < multiplicities[base]; ++m) 294 for (unsigned int local_index = 0; 295 local_index < fes[base]->n_dofs_per_quad(quad_number); 296 ++local_index, ++total_index) 297 { 298 const unsigned int index_in_base = 299 local_index + fes[base]->get_first_quad_index(quad_number); 300 301 Assert(index_in_base < fes[base]->n_dofs_per_cell(), 302 ExcInternalError()); 303 retval[total_index] = 304 fes[base]->restriction_is_additive(index_in_base); 305 } 306 } 307 308 // 4. Hexes 309 if (dim == 3) 310 { 311 for (unsigned int base = 0; base < fes.size(); ++base) 312 for (unsigned int m = 0; m < multiplicities[base]; ++m) 313 for (unsigned int local_index = 0; 314 local_index < fes[base]->n_dofs_per_hex(); 315 ++local_index, ++total_index) 316 { 317 const unsigned int index_in_base = 318 local_index + fes[base]->get_first_hex_index(); 319 320 Assert(index_in_base < fes[base]->n_dofs_per_cell(), 321 ExcInternalError()); 322 retval[total_index] = 323 fes[base]->restriction_is_additive(index_in_base); 324 } 325 } 326 327 Assert(total_index == n_shape_functions, ExcInternalError()); 328 329 return retval; 330 } 331 332 333 334 /** 335 * Take a @p FiniteElement object 336 * and return an boolean vector including the @p 337 * restriction_is_additive_flags of the mixed element consisting of @p N 338 * elements of the sub-element @p fe. 339 */ 340 template <int dim, int spacedim> 341 std::vector<bool> compute_restriction_is_additive_flags(const FiniteElement<dim,spacedim> * fe1,const unsigned int N1,const FiniteElement<dim,spacedim> * fe2,const unsigned int N2,const FiniteElement<dim,spacedim> * fe3,const unsigned int N3,const FiniteElement<dim,spacedim> * fe4,const unsigned int N4,const FiniteElement<dim,spacedim> * fe5,const unsigned int N5)342 compute_restriction_is_additive_flags( 343 const FiniteElement<dim, spacedim> *fe1, 344 const unsigned int N1, 345 const FiniteElement<dim, spacedim> *fe2, 346 const unsigned int N2, 347 const FiniteElement<dim, spacedim> *fe3, 348 const unsigned int N3, 349 const FiniteElement<dim, spacedim> *fe4, 350 const unsigned int N4, 351 const FiniteElement<dim, spacedim> *fe5, 352 const unsigned int N5) 353 { 354 std::vector<const FiniteElement<dim, spacedim> *> fe_list; 355 std::vector<unsigned int> multiplicities; 356 357 fe_list.push_back(fe1); 358 multiplicities.push_back(N1); 359 360 fe_list.push_back(fe2); 361 multiplicities.push_back(N2); 362 363 fe_list.push_back(fe3); 364 multiplicities.push_back(N3); 365 366 fe_list.push_back(fe4); 367 multiplicities.push_back(N4); 368 369 fe_list.push_back(fe5); 370 multiplicities.push_back(N5); 371 return compute_restriction_is_additive_flags(fe_list, multiplicities); 372 } 373 374 375 376 template <int dim, int spacedim> 377 std::vector<ComponentMask> compute_nonzero_components(const std::vector<const FiniteElement<dim,spacedim> * > & fes,const std::vector<unsigned int> & multiplicities,const bool do_tensor_product)378 compute_nonzero_components( 379 const std::vector<const FiniteElement<dim, spacedim> *> &fes, 380 const std::vector<unsigned int> & multiplicities, 381 const bool do_tensor_product) 382 { 383 AssertDimension(fes.size(), multiplicities.size()); 384 385 // first count the number of dofs and components that will emerge from the 386 // given FEs 387 unsigned int n_shape_functions = 0; 388 for (unsigned int i = 0; i < fes.size(); ++i) 389 if (multiplicities[i] > 0) // needed because fe might be nullptr 390 n_shape_functions += fes[i]->n_dofs_per_cell() * multiplicities[i]; 391 392 unsigned int n_components = 0; 393 if (do_tensor_product) 394 { 395 for (unsigned int i = 0; i < fes.size(); ++i) 396 if (multiplicities[i] > 0) // needed because fe might be nullptr 397 n_components += fes[i]->n_components() * multiplicities[i]; 398 } 399 else 400 { 401 for (unsigned int i = 0; i < fes.size(); ++i) 402 if (multiplicities[i] > 0) // needed because fe might be nullptr 403 { 404 n_components = fes[i]->n_components(); 405 break; 406 } 407 // Now check that all FEs have the same number of components: 408 for (unsigned int i = 0; i < fes.size(); ++i) 409 if (multiplicities[i] > 0) // needed because fe might be nullptr 410 Assert(n_components == fes[i]->n_components(), 411 ExcDimensionMismatch(n_components, 412 fes[i]->n_components())); 413 } 414 415 // generate the array that will hold the output 416 std::vector<std::vector<bool>> retval(n_shape_functions, 417 std::vector<bool>(n_components, 418 false)); 419 420 // finally go through all the shape functions of the base elements, and 421 // copy their flags. this somehow copies the code in build_cell_table, 422 // which is not nice as it uses too much implicit knowledge about the 423 // layout of the individual bases in the composed FE, but there seems no 424 // way around... 425 // 426 // for each shape function, copy the non-zero flags from the base element 427 // to this one, taking into account multiplicities, multiple components in 428 // base elements, and other complications 429 unsigned int total_index = 0; 430 for (const unsigned int vertex_number : 431 ReferenceCell::internal::Info::get_cell( 432 fes.front()->reference_cell_type()) 433 .vertex_indices()) 434 { 435 unsigned int comp_start = 0; 436 for (unsigned int base = 0; base < fes.size(); ++base) 437 for (unsigned int m = 0; m < multiplicities[base]; 438 ++m, 439 comp_start += 440 fes[base]->n_components() * do_tensor_product) 441 for (unsigned int local_index = 0; 442 local_index < fes[base]->n_dofs_per_vertex(); 443 ++local_index, ++total_index) 444 { 445 const unsigned int index_in_base = 446 (fes[base]->n_dofs_per_vertex() * vertex_number + 447 local_index); 448 449 Assert(comp_start + fes[base]->n_components() <= 450 retval[total_index].size(), 451 ExcInternalError()); 452 for (unsigned int c = 0; c < fes[base]->n_components(); ++c) 453 { 454 Assert(c < fes[base] 455 ->get_nonzero_components(index_in_base) 456 .size(), 457 ExcInternalError()); 458 retval[total_index][comp_start + c] = 459 fes[base]->get_nonzero_components(index_in_base)[c]; 460 } 461 } 462 } 463 464 // 2. Lines 465 for (const unsigned int line_number : 466 ReferenceCell::internal::Info::get_cell( 467 fes.front()->reference_cell_type()) 468 .line_indices()) 469 { 470 unsigned int comp_start = 0; 471 for (unsigned int base = 0; base < fes.size(); ++base) 472 for (unsigned int m = 0; m < multiplicities[base]; 473 ++m, 474 comp_start += 475 fes[base]->n_components() * do_tensor_product) 476 for (unsigned int local_index = 0; 477 local_index < fes[base]->n_dofs_per_line(); 478 ++local_index, ++total_index) 479 { 480 const unsigned int index_in_base = 481 (fes[base]->n_dofs_per_line() * line_number + local_index + 482 fes[base]->get_first_line_index()); 483 484 Assert(comp_start + fes[base]->n_components() <= 485 retval[total_index].size(), 486 ExcInternalError()); 487 for (unsigned int c = 0; c < fes[base]->n_components(); ++c) 488 { 489 Assert(c < fes[base] 490 ->get_nonzero_components(index_in_base) 491 .size(), 492 ExcInternalError()); 493 retval[total_index][comp_start + c] = 494 fes[base]->get_nonzero_components(index_in_base)[c]; 495 } 496 } 497 } 498 499 // 3. Quads 500 for (unsigned int quad_number = 0; 501 quad_number < (dim == 2 ? 502 1 : 503 (dim == 3 ? ReferenceCell::internal::Info::get_cell( 504 fes.front()->reference_cell_type()) 505 .n_faces() : 506 0)); 507 ++quad_number) 508 { 509 unsigned int comp_start = 0; 510 for (unsigned int base = 0; base < fes.size(); ++base) 511 for (unsigned int m = 0; m < multiplicities[base]; 512 ++m, 513 comp_start += 514 fes[base]->n_components() * do_tensor_product) 515 for (unsigned int local_index = 0; 516 local_index < fes[base]->n_dofs_per_quad(quad_number); 517 ++local_index, ++total_index) 518 { 519 const unsigned int index_in_base = 520 local_index + fes[base]->get_first_quad_index(quad_number); 521 522 Assert(comp_start + fes[base]->n_components() <= 523 retval[total_index].size(), 524 ExcInternalError()); 525 for (unsigned int c = 0; c < fes[base]->n_components(); ++c) 526 { 527 Assert(c < fes[base] 528 ->get_nonzero_components(index_in_base) 529 .size(), 530 ExcInternalError()); 531 retval[total_index][comp_start + c] = 532 fes[base]->get_nonzero_components(index_in_base)[c]; 533 } 534 } 535 } 536 537 // 4. Hexes 538 if (dim == 3) 539 { 540 unsigned int comp_start = 0; 541 for (unsigned int base = 0; base < fes.size(); ++base) 542 for (unsigned int m = 0; m < multiplicities[base]; 543 ++m, 544 comp_start += 545 fes[base]->n_components() * do_tensor_product) 546 for (unsigned int local_index = 0; 547 local_index < fes[base]->n_dofs_per_hex(); 548 ++local_index, ++total_index) 549 { 550 const unsigned int index_in_base = 551 local_index + fes[base]->get_first_hex_index(); 552 553 Assert(comp_start + fes[base]->n_components() <= 554 retval[total_index].size(), 555 ExcInternalError()); 556 for (unsigned int c = 0; c < fes[base]->n_components(); ++c) 557 { 558 Assert(c < fes[base] 559 ->get_nonzero_components(index_in_base) 560 .size(), 561 ExcInternalError()); 562 retval[total_index][comp_start + c] = 563 fes[base]->get_nonzero_components(index_in_base)[c]; 564 } 565 } 566 } 567 568 Assert(total_index == n_shape_functions, ExcInternalError()); 569 570 // now copy the vector<vector<bool> > into a vector<ComponentMask>. 571 // this appears complicated but we do it this way since it's just 572 // awkward to generate ComponentMasks directly and so we need the 573 // recourse of the inner vector<bool> anyway. 574 std::vector<ComponentMask> xretval(retval.size()); 575 for (unsigned int i = 0; i < retval.size(); ++i) 576 xretval[i] = ComponentMask(retval[i]); 577 return xretval; 578 } 579 580 581 582 /** 583 * Compute the non-zero vector components of a composed finite element. 584 */ 585 template <int dim, int spacedim> 586 std::vector<ComponentMask> compute_nonzero_components(const FiniteElement<dim,spacedim> * fe1,const unsigned int N1,const FiniteElement<dim,spacedim> * fe2,const unsigned int N2,const FiniteElement<dim,spacedim> * fe3,const unsigned int N3,const FiniteElement<dim,spacedim> * fe4,const unsigned int N4,const FiniteElement<dim,spacedim> * fe5,const unsigned int N5,const bool do_tensor_product)587 compute_nonzero_components(const FiniteElement<dim, spacedim> *fe1, 588 const unsigned int N1, 589 const FiniteElement<dim, spacedim> *fe2, 590 const unsigned int N2, 591 const FiniteElement<dim, spacedim> *fe3, 592 const unsigned int N3, 593 const FiniteElement<dim, spacedim> *fe4, 594 const unsigned int N4, 595 const FiniteElement<dim, spacedim> *fe5, 596 const unsigned int N5, 597 const bool do_tensor_product) 598 { 599 std::vector<const FiniteElement<dim, spacedim> *> fe_list; 600 std::vector<unsigned int> multiplicities; 601 602 fe_list.push_back(fe1); 603 multiplicities.push_back(N1); 604 605 fe_list.push_back(fe2); 606 multiplicities.push_back(N2); 607 608 fe_list.push_back(fe3); 609 multiplicities.push_back(N3); 610 611 fe_list.push_back(fe4); 612 multiplicities.push_back(N4); 613 614 fe_list.push_back(fe5); 615 multiplicities.push_back(N5); 616 617 return compute_nonzero_components(fe_list, 618 multiplicities, 619 do_tensor_product); 620 } 621 622 623 624 template <int dim, int spacedim> 625 void build_cell_tables(std::vector<std::pair<std::pair<unsigned int,unsigned int>,unsigned int>> & system_to_base_table,std::vector<std::pair<unsigned int,unsigned int>> & system_to_component_table,std::vector<std::pair<std::pair<unsigned int,unsigned int>,unsigned int>> & component_to_base_table,const FiniteElement<dim,spacedim> & fe,const bool do_tensor_product)626 build_cell_tables( 627 std::vector<std::pair<std::pair<unsigned int, unsigned int>, 628 unsigned int>> &system_to_base_table, 629 std::vector<std::pair<unsigned int, unsigned int>> 630 & system_to_component_table, 631 std::vector<std::pair<std::pair<unsigned int, unsigned int>, 632 unsigned int>> &component_to_base_table, 633 const FiniteElement<dim, spacedim> & fe, 634 const bool do_tensor_product) 635 { 636 unsigned int total_index = 0; 637 638 if (do_tensor_product) 639 { 640 for (unsigned int base = 0; base < fe.n_base_elements(); ++base) 641 for (unsigned int m = 0; m < fe.element_multiplicity(base); ++m) 642 { 643 for (unsigned int k = 0; 644 k < fe.base_element(base).n_components(); 645 ++k) 646 component_to_base_table[total_index++] = 647 std::make_pair(std::make_pair(base, k), m); 648 } 649 Assert(total_index == component_to_base_table.size(), 650 ExcInternalError()); 651 } 652 else 653 { 654 // The base element establishing a component does not make sense in 655 // this case. Set up to something meaningless: 656 std::fill( 657 component_to_base_table.begin(), 658 component_to_base_table.end(), 659 std::make_pair(std::make_pair(numbers::invalid_unsigned_int, 660 numbers::invalid_unsigned_int), 661 numbers::invalid_unsigned_int)); 662 } 663 664 665 // Initialize index tables. Multi-component base elements have to be 666 // thought of. For non-primitive shape functions, have a special invalid 667 // index. 668 const std::pair<unsigned int, unsigned int> non_primitive_index( 669 numbers::invalid_unsigned_int, numbers::invalid_unsigned_int); 670 671 // First enumerate vertex indices, where we first enumerate all indices on 672 // the first vertex in the order of the base elements, then of the second 673 // vertex, etc 674 total_index = 0; 675 for (const unsigned int vertex_number : 676 ReferenceCell::internal::Info::get_cell(fe.reference_cell_type()) 677 .vertex_indices()) 678 { 679 unsigned int comp_start = 0; 680 for (unsigned int base = 0; base < fe.n_base_elements(); ++base) 681 for (unsigned int m = 0; m < fe.element_multiplicity(base); 682 ++m, 683 comp_start += 684 fe.base_element(base).n_components() * 685 do_tensor_product) 686 for (unsigned int local_index = 0; 687 local_index < fe.base_element(base).n_dofs_per_vertex(); 688 ++local_index, ++total_index) 689 { 690 const unsigned int index_in_base = 691 (fe.base_element(base).n_dofs_per_vertex() * vertex_number + 692 local_index); 693 694 system_to_base_table[total_index] = 695 std::make_pair(std::make_pair(base, m), index_in_base); 696 697 if (fe.base_element(base).is_primitive(index_in_base)) 698 { 699 const unsigned int comp_in_base = 700 fe.base_element(base) 701 .system_to_component_index(index_in_base) 702 .first; 703 const unsigned int comp = comp_start + comp_in_base; 704 const unsigned int index_in_comp = 705 fe.base_element(base) 706 .system_to_component_index(index_in_base) 707 .second; 708 system_to_component_table[total_index] = 709 std::make_pair(comp, index_in_comp); 710 } 711 else 712 system_to_component_table[total_index] = 713 non_primitive_index; 714 } 715 } 716 717 // 2. Lines 718 for (const unsigned int line_number : 719 ReferenceCell::internal::Info::get_cell(fe.reference_cell_type()) 720 .line_indices()) 721 { 722 unsigned int comp_start = 0; 723 for (unsigned int base = 0; base < fe.n_base_elements(); ++base) 724 for (unsigned int m = 0; m < fe.element_multiplicity(base); 725 ++m, 726 comp_start += 727 fe.base_element(base).n_components() * 728 do_tensor_product) 729 for (unsigned int local_index = 0; 730 local_index < fe.base_element(base).n_dofs_per_line(); 731 ++local_index, ++total_index) 732 { 733 const unsigned int index_in_base = 734 (fe.base_element(base).n_dofs_per_line() * line_number + 735 local_index + 736 fe.base_element(base).get_first_line_index()); 737 738 system_to_base_table[total_index] = 739 std::make_pair(std::make_pair(base, m), index_in_base); 740 741 if (fe.base_element(base).is_primitive(index_in_base)) 742 { 743 const unsigned int comp_in_base = 744 fe.base_element(base) 745 .system_to_component_index(index_in_base) 746 .first; 747 const unsigned int comp = comp_start + comp_in_base; 748 const unsigned int index_in_comp = 749 fe.base_element(base) 750 .system_to_component_index(index_in_base) 751 .second; 752 system_to_component_table[total_index] = 753 std::make_pair(comp, index_in_comp); 754 } 755 else 756 system_to_component_table[total_index] = 757 non_primitive_index; 758 } 759 } 760 761 // 3. Quads 762 for (unsigned int quad_number = 0; 763 quad_number < (dim == 2 ? 764 1 : 765 (dim == 3 ? ReferenceCell::internal::Info::get_cell( 766 fe.reference_cell_type()) 767 .n_faces() : 768 0)); 769 ++quad_number) 770 { 771 unsigned int comp_start = 0; 772 for (unsigned int base = 0; base < fe.n_base_elements(); ++base) 773 for (unsigned int m = 0; m < fe.element_multiplicity(base); 774 ++m, 775 comp_start += 776 fe.base_element(base).n_components() * 777 do_tensor_product) 778 for (unsigned int local_index = 0; 779 local_index < 780 fe.base_element(base).n_dofs_per_quad(quad_number); 781 ++local_index, ++total_index) 782 { 783 const unsigned int index_in_base = 784 local_index + 785 fe.base_element(base).get_first_quad_index(quad_number); 786 787 system_to_base_table[total_index] = 788 std::make_pair(std::make_pair(base, m), index_in_base); 789 790 if (fe.base_element(base).is_primitive(index_in_base)) 791 { 792 const unsigned int comp_in_base = 793 fe.base_element(base) 794 .system_to_component_index(index_in_base) 795 .first; 796 const unsigned int comp = comp_start + comp_in_base; 797 const unsigned int index_in_comp = 798 fe.base_element(base) 799 .system_to_component_index(index_in_base) 800 .second; 801 system_to_component_table[total_index] = 802 std::make_pair(comp, index_in_comp); 803 } 804 else 805 system_to_component_table[total_index] = 806 non_primitive_index; 807 } 808 } 809 810 // 4. Hexes 811 if (dim == 3) 812 { 813 unsigned int comp_start = 0; 814 for (unsigned int base = 0; base < fe.n_base_elements(); ++base) 815 for (unsigned int m = 0; m < fe.element_multiplicity(base); 816 ++m, 817 comp_start += 818 fe.base_element(base).n_components() * 819 do_tensor_product) 820 for (unsigned int local_index = 0; 821 local_index < fe.base_element(base).n_dofs_per_hex(); 822 ++local_index, ++total_index) 823 { 824 const unsigned int index_in_base = 825 local_index + fe.base_element(base).get_first_hex_index(); 826 827 system_to_base_table[total_index] = 828 std::make_pair(std::make_pair(base, m), index_in_base); 829 830 if (fe.base_element(base).is_primitive(index_in_base)) 831 { 832 const unsigned int comp_in_base = 833 fe.base_element(base) 834 .system_to_component_index(index_in_base) 835 .first; 836 const unsigned int comp = comp_start + comp_in_base; 837 const unsigned int index_in_comp = 838 fe.base_element(base) 839 .system_to_component_index(index_in_base) 840 .second; 841 system_to_component_table[total_index] = 842 std::make_pair(comp, index_in_comp); 843 } 844 else 845 system_to_component_table[total_index] = 846 non_primitive_index; 847 } 848 } 849 } 850 851 852 853 template <int dim, int spacedim> 854 void build_face_tables(std::vector<std::pair<std::pair<unsigned int,unsigned int>,unsigned int>> & face_system_to_base_table,std::vector<std::pair<unsigned int,unsigned int>> & face_system_to_component_table,const FiniteElement<dim,spacedim> & fe,const bool do_tensor_product,const unsigned int face_no)855 build_face_tables( 856 std::vector<std::pair<std::pair<unsigned int, unsigned int>, 857 unsigned int>> &face_system_to_base_table, 858 std::vector<std::pair<unsigned int, unsigned int>> 859 & face_system_to_component_table, 860 const FiniteElement<dim, spacedim> &fe, 861 const bool do_tensor_product, 862 const unsigned int face_no) 863 { 864 // Initialize index tables. do this in the same way as done for the cell 865 // tables, except that we now loop over the objects of faces 866 867 // For non-primitive shape functions, have a special invalid index 868 const std::pair<unsigned int, unsigned int> non_primitive_index( 869 numbers::invalid_unsigned_int, numbers::invalid_unsigned_int); 870 871 // 1. Vertices 872 unsigned int total_index = 0; 873 for (unsigned int vertex_number = 0; 874 vertex_number < 875 ReferenceCell::internal::Info::get_face(fe.reference_cell_type(), 876 face_no) 877 .n_vertices(); 878 ++vertex_number) 879 { 880 unsigned int comp_start = 0; 881 for (unsigned int base = 0; base < fe.n_base_elements(); ++base) 882 for (unsigned int m = 0; m < fe.element_multiplicity(base); 883 ++m, 884 comp_start += 885 fe.base_element(base).n_components() * 886 do_tensor_product) 887 for (unsigned int local_index = 0; 888 local_index < fe.base_element(base).n_dofs_per_vertex(); 889 ++local_index, ++total_index) 890 { 891 // get (cell) index of this shape function inside the base 892 // element to see whether the shape function is primitive 893 // (assume that all shape functions on vertices share the same 894 // primitivity property; assume likewise for all shape 895 // functions located on lines, quads, etc. this way, we can 896 // ask for primitivity of only _one_ shape function, which is 897 // taken as representative for all others located on the same 898 // type of object): 899 const unsigned int index_in_base = 900 (fe.base_element(base).n_dofs_per_vertex() * vertex_number + 901 local_index); 902 903 const unsigned int face_index_in_base = 904 (fe.base_element(base).n_dofs_per_vertex() * vertex_number + 905 local_index); 906 907 face_system_to_base_table[total_index] = 908 std::make_pair(std::make_pair(base, m), face_index_in_base); 909 910 if (fe.base_element(base).is_primitive(index_in_base)) 911 { 912 const unsigned int comp_in_base = 913 fe.base_element(base) 914 .face_system_to_component_index(face_index_in_base, 915 face_no) 916 .first; 917 const unsigned int comp = comp_start + comp_in_base; 918 const unsigned int face_index_in_comp = 919 fe.base_element(base) 920 .face_system_to_component_index(face_index_in_base, 921 face_no) 922 .second; 923 face_system_to_component_table[total_index] = 924 std::make_pair(comp, face_index_in_comp); 925 } 926 else 927 face_system_to_component_table[total_index] = 928 non_primitive_index; 929 } 930 } 931 932 // 2. Lines 933 for (unsigned int line_number = 0; 934 line_number < 935 ReferenceCell::internal::Info::get_face(fe.reference_cell_type(), 936 face_no) 937 .n_lines(); 938 ++line_number) 939 { 940 unsigned int comp_start = 0; 941 for (unsigned int base = 0; base < fe.n_base_elements(); ++base) 942 for (unsigned int m = 0; m < fe.element_multiplicity(base); 943 ++m, 944 comp_start += 945 fe.base_element(base).n_components() * 946 do_tensor_product) 947 for (unsigned int local_index = 0; 948 local_index < fe.base_element(base).n_dofs_per_line(); 949 ++local_index, ++total_index) 950 { 951 // do everything alike for this type of object 952 const unsigned int index_in_base = 953 (fe.base_element(base).n_dofs_per_line() * line_number + 954 local_index + 955 fe.base_element(base).get_first_line_index()); 956 957 const unsigned int face_index_in_base = 958 (fe.base_element(base).get_first_face_line_index(face_no) + 959 fe.base_element(base).n_dofs_per_line() * line_number + 960 local_index); 961 962 face_system_to_base_table[total_index] = 963 std::make_pair(std::make_pair(base, m), face_index_in_base); 964 965 if (fe.base_element(base).is_primitive(index_in_base)) 966 { 967 const unsigned int comp_in_base = 968 fe.base_element(base) 969 .face_system_to_component_index(face_index_in_base, 970 face_no) 971 .first; 972 const unsigned int comp = comp_start + comp_in_base; 973 const unsigned int face_index_in_comp = 974 fe.base_element(base) 975 .face_system_to_component_index(face_index_in_base, 976 face_no) 977 .second; 978 face_system_to_component_table[total_index] = 979 std::make_pair(comp, face_index_in_comp); 980 } 981 else 982 face_system_to_component_table[total_index] = 983 non_primitive_index; 984 } 985 } 986 987 // 3. Quads 988 if (dim == 3) 989 { 990 unsigned int comp_start = 0; 991 for (unsigned int base = 0; base < fe.n_base_elements(); ++base) 992 for (unsigned int m = 0; m < fe.element_multiplicity(base); 993 ++m, 994 comp_start += 995 fe.base_element(base).n_components() * 996 do_tensor_product) 997 for (unsigned int local_index = 0; 998 local_index < fe.base_element(base).n_dofs_per_quad(face_no); 999 ++local_index, ++total_index) 1000 { 1001 // do everything alike for this type of object 1002 const unsigned int index_in_base = 1003 (local_index + 1004 fe.base_element(base).get_first_quad_index(face_no)); 1005 1006 const unsigned int face_index_in_base = 1007 (fe.base_element(base).get_first_face_quad_index(face_no) + 1008 local_index); 1009 1010 face_system_to_base_table[total_index] = 1011 std::make_pair(std::make_pair(base, m), face_index_in_base); 1012 1013 if (fe.base_element(base).is_primitive(index_in_base)) 1014 { 1015 const unsigned int comp_in_base = 1016 fe.base_element(base) 1017 .face_system_to_component_index(face_index_in_base, 1018 face_no) 1019 .first; 1020 const unsigned int comp = comp_start + comp_in_base; 1021 const unsigned int face_index_in_comp = 1022 fe.base_element(base) 1023 .face_system_to_component_index(face_index_in_base, 1024 face_no) 1025 .second; 1026 face_system_to_component_table[total_index] = 1027 std::make_pair(comp, face_index_in_comp); 1028 } 1029 else 1030 face_system_to_component_table[total_index] = 1031 non_primitive_index; 1032 } 1033 } 1034 Assert(total_index == fe.n_dofs_per_face(face_no), ExcInternalError()); 1035 Assert(total_index == face_system_to_component_table.size(), 1036 ExcInternalError()); 1037 Assert(total_index == face_system_to_base_table.size(), 1038 ExcInternalError()); 1039 } 1040 } // namespace Compositing 1041 1042 1043 1044 // Not implemented in the general case. 1045 template <class FE> 1046 std::unique_ptr<FiniteElement<FE::dimension, FE::space_dimension>> get(const Quadrature<1> &)1047 FEFactory<FE>::get(const Quadrature<1> &) const 1048 { 1049 Assert(false, ExcNotImplemented()); 1050 return nullptr; 1051 } 1052 1053 // Specializations for FE_Q. 1054 template <> 1055 std::unique_ptr<FiniteElement<1, 1>> get(const Quadrature<1> & quad)1056 FEFactory<FE_Q<1, 1>>::get(const Quadrature<1> &quad) const 1057 { 1058 return std::make_unique<FE_Q<1>>(quad); 1059 } 1060 1061 template <> 1062 std::unique_ptr<FiniteElement<2, 2>> get(const Quadrature<1> & quad)1063 FEFactory<FE_Q<2, 2>>::get(const Quadrature<1> &quad) const 1064 { 1065 return std::make_unique<FE_Q<2>>(quad); 1066 } 1067 1068 template <> 1069 std::unique_ptr<FiniteElement<3, 3>> get(const Quadrature<1> & quad)1070 FEFactory<FE_Q<3, 3>>::get(const Quadrature<1> &quad) const 1071 { 1072 return std::make_unique<FE_Q<3>>(quad); 1073 } 1074 1075 // Specializations for FE_Q_DG0. 1076 template <> 1077 std::unique_ptr<FiniteElement<1, 1>> get(const Quadrature<1> & quad)1078 FEFactory<FE_Q_DG0<1, 1>>::get(const Quadrature<1> &quad) const 1079 { 1080 return std::make_unique<FE_Q_DG0<1>>(quad); 1081 } 1082 1083 template <> 1084 std::unique_ptr<FiniteElement<2, 2>> get(const Quadrature<1> & quad)1085 FEFactory<FE_Q_DG0<2, 2>>::get(const Quadrature<1> &quad) const 1086 { 1087 return std::make_unique<FE_Q_DG0<2>>(quad); 1088 } 1089 1090 template <> 1091 std::unique_ptr<FiniteElement<3, 3>> get(const Quadrature<1> & quad)1092 FEFactory<FE_Q_DG0<3, 3>>::get(const Quadrature<1> &quad) const 1093 { 1094 return std::make_unique<FE_Q_DG0<3>>(quad); 1095 } 1096 1097 // Specializations for FE_Q_Bubbles. 1098 template <> 1099 std::unique_ptr<FiniteElement<1, 1>> get(const Quadrature<1> & quad)1100 FEFactory<FE_Q_Bubbles<1, 1>>::get(const Quadrature<1> &quad) const 1101 { 1102 return std::make_unique<FE_Q_Bubbles<1>>(quad); 1103 } 1104 1105 template <> 1106 std::unique_ptr<FiniteElement<2, 2>> get(const Quadrature<1> & quad)1107 FEFactory<FE_Q_Bubbles<2, 2>>::get(const Quadrature<1> &quad) const 1108 { 1109 return std::make_unique<FE_Q_Bubbles<2>>(quad); 1110 } 1111 1112 template <> 1113 std::unique_ptr<FiniteElement<3, 3>> get(const Quadrature<1> & quad)1114 FEFactory<FE_Q_Bubbles<3, 3>>::get(const Quadrature<1> &quad) const 1115 { 1116 return std::make_unique<FE_Q_Bubbles<3>>(quad); 1117 } 1118 1119 // Specializations for FE_DGQArbitraryNodes. 1120 template <> 1121 std::unique_ptr<FiniteElement<1, 1>> get(const Quadrature<1> & quad)1122 FEFactory<FE_DGQ<1>>::get(const Quadrature<1> &quad) const 1123 { 1124 return std::make_unique<FE_DGQArbitraryNodes<1>>(quad); 1125 } 1126 1127 template <> 1128 std::unique_ptr<FiniteElement<1, 2>> get(const Quadrature<1> & quad)1129 FEFactory<FE_DGQ<1, 2>>::get(const Quadrature<1> &quad) const 1130 { 1131 return std::make_unique<FE_DGQArbitraryNodes<1, 2>>(quad); 1132 } 1133 1134 template <> 1135 std::unique_ptr<FiniteElement<1, 3>> get(const Quadrature<1> & quad)1136 FEFactory<FE_DGQ<1, 3>>::get(const Quadrature<1> &quad) const 1137 { 1138 return std::make_unique<FE_DGQArbitraryNodes<1, 3>>(quad); 1139 } 1140 1141 template <> 1142 std::unique_ptr<FiniteElement<2, 2>> get(const Quadrature<1> & quad)1143 FEFactory<FE_DGQ<2>>::get(const Quadrature<1> &quad) const 1144 { 1145 return std::make_unique<FE_DGQArbitraryNodes<2>>(quad); 1146 } 1147 1148 template <> 1149 std::unique_ptr<FiniteElement<2, 3>> get(const Quadrature<1> & quad)1150 FEFactory<FE_DGQ<2, 3>>::get(const Quadrature<1> &quad) const 1151 { 1152 return std::make_unique<FE_DGQArbitraryNodes<2, 3>>(quad); 1153 } 1154 1155 template <> 1156 std::unique_ptr<FiniteElement<3, 3>> get(const Quadrature<1> & quad)1157 FEFactory<FE_DGQ<3>>::get(const Quadrature<1> &quad) const 1158 { 1159 return std::make_unique<FE_DGQArbitraryNodes<3>>(quad); 1160 } 1161 1162 1163 1164 namespace internal 1165 { 1166 // The following three functions serve to fill the maps from element 1167 // names to elements fe_name_map below. The first one exists because 1168 // we have finite elements which are not implemented for nonzero 1169 // codimension. These should be transferred to the second function 1170 // eventually. 1171 namespace FEToolsAddFENameHelper 1172 { 1173 template <int dim> 1174 void fill_no_codim_fe_names(std::map<std::string,std::unique_ptr<const Subscriptor>> & result)1175 fill_no_codim_fe_names( 1176 std::map<std::string, std::unique_ptr<const Subscriptor>> &result) 1177 { 1178 result["FE_Q_Hierarchical"] = 1179 std::make_unique<FETools::FEFactory<FE_Q_Hierarchical<dim>>>(); 1180 result["FE_ABF"] = std::make_unique<FETools::FEFactory<FE_ABF<dim>>>(); 1181 result["FE_Bernstein"] = 1182 std::make_unique<FETools::FEFactory<FE_Bernstein<dim>>>(); 1183 result["FE_BDM"] = std::make_unique<FETools::FEFactory<FE_BDM<dim>>>(); 1184 result["FE_DGBDM"] = 1185 std::make_unique<FETools::FEFactory<FE_DGBDM<dim>>>(); 1186 result["FE_DGNedelec"] = 1187 std::make_unique<FETools::FEFactory<FE_DGNedelec<dim>>>(); 1188 result["FE_DGRaviartThomas"] = 1189 std::make_unique<FETools::FEFactory<FE_DGRaviartThomas<dim>>>(); 1190 result["FE_RaviartThomas"] = 1191 std::make_unique<FETools::FEFactory<FE_RaviartThomas<dim>>>(); 1192 result["FE_RaviartThomasNodal"] = 1193 std::make_unique<FETools::FEFactory<FE_RaviartThomasNodal<dim>>>(); 1194 result["FE_RT_Bubbles"] = 1195 std::make_unique<FETools::FEFactory<FE_RT_Bubbles<dim>>>(); 1196 result["FE_Nedelec"] = 1197 std::make_unique<FETools::FEFactory<FE_Nedelec<dim>>>(); 1198 result["FE_NedelecSZ"] = 1199 std::make_unique<FETools::FEFactory<FE_NedelecSZ<dim>>>(); 1200 result["FE_DGPNonparametric"] = 1201 std::make_unique<FETools::FEFactory<FE_DGPNonparametric<dim>>>(); 1202 result["FE_DGP"] = std::make_unique<FETools::FEFactory<FE_DGP<dim>>>(); 1203 result["FE_DGPMonomial"] = 1204 std::make_unique<FETools::FEFactory<FE_DGPMonomial<dim>>>(); 1205 result["FE_DGQ"] = std::make_unique<FETools::FEFactory<FE_DGQ<dim>>>(); 1206 result["FE_DGQArbitraryNodes"] = 1207 std::make_unique<FETools::FEFactory<FE_DGQ<dim>>>(); 1208 result["FE_DGQLegendre"] = 1209 std::make_unique<FETools::FEFactory<FE_DGQLegendre<dim>>>(); 1210 result["FE_DGQHermite"] = 1211 std::make_unique<FETools::FEFactory<FE_DGQHermite<dim>>>(); 1212 result["FE_FaceQ"] = 1213 std::make_unique<FETools::FEFactory<FE_FaceQ<dim>>>(); 1214 result["FE_FaceP"] = 1215 std::make_unique<FETools::FEFactory<FE_FaceP<dim>>>(); 1216 result["FE_Q"] = std::make_unique<FETools::FEFactory<FE_Q<dim>>>(); 1217 result["FE_Q_DG0"] = 1218 std::make_unique<FETools::FEFactory<FE_Q_DG0<dim>>>(); 1219 result["FE_Q_Bubbles"] = 1220 std::make_unique<FETools::FEFactory<FE_Q_Bubbles<dim>>>(); 1221 result["FE_Q_iso_Q1"] = 1222 std::make_unique<FETools::FEFactory<FE_Q_iso_Q1<dim>>>(); 1223 result["FE_Nothing"] = 1224 std::make_unique<FETools::FEFactory<FE_Nothing<dim>>>(); 1225 result["FE_RannacherTurek"] = 1226 std::make_unique<FETools::FEFactory<FE_RannacherTurek<dim>>>(); 1227 } 1228 1229 1230 1231 // This function fills a map from names to finite elements for any 1232 // dimension and codimension for those elements which support 1233 // nonzero codimension. 1234 template <int dim, int spacedim> 1235 void fill_codim_fe_names(std::map<std::string,std::unique_ptr<const Subscriptor>> & result)1236 fill_codim_fe_names( 1237 std::map<std::string, std::unique_ptr<const Subscriptor>> &result) 1238 { 1239 result["FE_Bernstein"] = 1240 std::make_unique<FETools::FEFactory<FE_Bernstein<dim, spacedim>>>(); 1241 result["FE_DGP"] = 1242 std::make_unique<FETools::FEFactory<FE_DGP<dim, spacedim>>>(); 1243 result["FE_DGQ"] = 1244 std::make_unique<FETools::FEFactory<FE_DGQ<dim, spacedim>>>(); 1245 result["FE_Nothing"] = 1246 std::make_unique<FETools::FEFactory<FE_Nothing<dim, spacedim>>>(); 1247 result["FE_DGQArbitraryNodes"] = 1248 std::make_unique<FETools::FEFactory<FE_DGQ<dim, spacedim>>>(); 1249 result["FE_DGQLegendre"] = 1250 std::make_unique<FETools::FEFactory<FE_DGQLegendre<dim, spacedim>>>(); 1251 result["FE_DGQHermite"] = 1252 std::make_unique<FETools::FEFactory<FE_DGQHermite<dim, spacedim>>>(); 1253 result["FE_Q_Bubbles"] = 1254 std::make_unique<FETools::FEFactory<FE_Q_Bubbles<dim, spacedim>>>(); 1255 result["FE_Q_DG0"] = 1256 std::make_unique<FETools::FEFactory<FE_Q_DG0<dim, spacedim>>>(); 1257 result["FE_Q_iso_Q1"] = 1258 std::make_unique<FETools::FEFactory<FE_Q_iso_Q1<dim, spacedim>>>(); 1259 result["FE_Q"] = 1260 std::make_unique<FETools::FEFactory<FE_Q<dim, spacedim>>>(); 1261 result["FE_Bernstein"] = 1262 std::make_unique<FETools::FEFactory<FE_Bernstein<dim, spacedim>>>(); 1263 } 1264 1265 // The function filling the vector fe_name_map below. It iterates 1266 // through all legal dimension/spacedimension pairs and fills 1267 // fe_name_map[dimension][spacedimension] with the maps generated 1268 // by the functions above. 1269 std::array< 1270 std::array<std::map<std::string, std::unique_ptr<const Subscriptor>>, 1271 4>, 1272 4> fill_default_map()1273 fill_default_map() 1274 { 1275 std::array< 1276 std::array<std::map<std::string, std::unique_ptr<const Subscriptor>>, 1277 4>, 1278 4> 1279 result; 1280 1281 fill_no_codim_fe_names<1>(result[1][1]); 1282 fill_no_codim_fe_names<2>(result[2][2]); 1283 fill_no_codim_fe_names<3>(result[3][3]); 1284 1285 fill_codim_fe_names<1, 2>(result[1][2]); 1286 fill_codim_fe_names<1, 3>(result[1][3]); 1287 fill_codim_fe_names<2, 3>(result[2][3]); 1288 1289 return result; 1290 } 1291 1292 1293 // have a lock that guarantees that at most one thread is changing 1294 // and accessing the fe_name_map variable. make this lock local to 1295 // this file. 1296 // 1297 // This variable is declared static (even though 1298 // it belongs to an internal namespace) in order to make icc happy 1299 // (which otherwise reports a multiply defined symbol when linking 1300 // libraries for more than one space dimension together) 1301 static std::mutex fe_name_map_lock; 1302 1303 // This is the map used by FETools::get_fe_by_name and 1304 // FETools::add_fe_name. It is only accessed by functions in this 1305 // file, so it is safe to make it a static variable here. It must be 1306 // static so that we can link several dimensions together. 1307 1308 // The organization of this storage is such that 1309 // fe_name_map[dim][spacedim][name] points to an 1310 // FEFactoryBase<dim,spacedim> with the name given. Since 1311 // all entries of this vector are of different type, we store 1312 // pointers to generic objects and cast them when needed. 1313 1314 // We use a unique pointer to factory objects, to ensure that they 1315 // get deleted at the end of the program run and don't end up as 1316 // apparent memory leaks to programs like valgrind. 1317 1318 // This vector is initialized at program start time using the 1319 // function above. because at this time there are no threads 1320 // running, there are no thread-safety issues here. since this is 1321 // compiled for all dimensions at once, need to create objects for 1322 // each dimension and then separate between them further down 1323 inline std::array< 1324 std::array<std::map<std::string, std::unique_ptr<const Subscriptor>>, 1325 4>, 1326 4> & get_fe_name_map()1327 get_fe_name_map() 1328 { 1329 static std::array< 1330 std::array<std::map<std::string, std::unique_ptr<const Subscriptor>>, 1331 4>, 1332 4> 1333 fe_name_map = fill_default_map(); 1334 return fe_name_map; 1335 } 1336 } // namespace FEToolsAddFENameHelper 1337 1338 namespace FEToolsGetInterpolationMatrixHelper 1339 { 1340 // forwarder function for 1341 // FE::get_interpolation_matrix. we 1342 // will want to call that function 1343 // for arbitrary FullMatrix<T> 1344 // types, but it only accepts 1345 // double arguments. since it is a 1346 // virtual function, this can also 1347 // not be changed. so have a 1348 // forwarder function that calls 1349 // that function directly if 1350 // T==double, and otherwise uses a 1351 // temporary 1352 template <int dim, int spacedim> 1353 inline void gim_forwarder(const FiniteElement<dim,spacedim> & fe1,const FiniteElement<dim,spacedim> & fe2,FullMatrix<double> & interpolation_matrix)1354 gim_forwarder(const FiniteElement<dim, spacedim> &fe1, 1355 const FiniteElement<dim, spacedim> &fe2, 1356 FullMatrix<double> & interpolation_matrix) 1357 { 1358 fe2.get_interpolation_matrix(fe1, interpolation_matrix); 1359 } 1360 1361 1362 1363 template <int dim, typename number, int spacedim> 1364 inline void gim_forwarder(const FiniteElement<dim,spacedim> & fe1,const FiniteElement<dim,spacedim> & fe2,FullMatrix<number> & interpolation_matrix)1365 gim_forwarder(const FiniteElement<dim, spacedim> &fe1, 1366 const FiniteElement<dim, spacedim> &fe2, 1367 FullMatrix<number> & interpolation_matrix) 1368 { 1369 FullMatrix<double> tmp(interpolation_matrix.m(), 1370 interpolation_matrix.n()); 1371 fe2.get_interpolation_matrix(fe1, tmp); 1372 interpolation_matrix = tmp; 1373 } 1374 } // namespace FEToolsGetInterpolationMatrixHelper 1375 } // namespace internal 1376 1377 1378 1379 template <int dim, int spacedim> 1380 void compute_component_wise(const FiniteElement<dim,spacedim> & element,std::vector<unsigned int> & renumbering,std::vector<std::vector<unsigned int>> & comp_start)1381 compute_component_wise(const FiniteElement<dim, spacedim> & element, 1382 std::vector<unsigned int> & renumbering, 1383 std::vector<std::vector<unsigned int>> &comp_start) 1384 { 1385 Assert(renumbering.size() == element.n_dofs_per_cell(), 1386 ExcDimensionMismatch(renumbering.size(), element.n_dofs_per_cell())); 1387 1388 comp_start.resize(element.n_base_elements()); 1389 1390 unsigned int index = 0; 1391 for (unsigned int i = 0; i < comp_start.size(); ++i) 1392 { 1393 comp_start[i].resize(element.element_multiplicity(i)); 1394 const unsigned int increment = 1395 element.base_element(i).n_dofs_per_cell(); 1396 1397 for (unsigned int &first_index_of_component : comp_start[i]) 1398 { 1399 first_index_of_component = index; 1400 index += increment; 1401 } 1402 } 1403 1404 // For each index i of the 1405 // unstructured cellwise 1406 // numbering, renumbering 1407 // contains the index of the 1408 // cell-block numbering 1409 for (unsigned int i = 0; i < element.n_dofs_per_cell(); ++i) 1410 { 1411 std::pair<std::pair<unsigned int, unsigned int>, unsigned int> indices = 1412 element.system_to_base_index(i); 1413 renumbering[i] = comp_start[indices.first.first][indices.first.second] + 1414 indices.second; 1415 } 1416 } 1417 1418 1419 1420 template <int dim, int spacedim> 1421 void compute_block_renumbering(const FiniteElement<dim,spacedim> & element,std::vector<types::global_dof_index> & renumbering,std::vector<types::global_dof_index> & block_data,bool return_start_indices)1422 compute_block_renumbering(const FiniteElement<dim, spacedim> & element, 1423 std::vector<types::global_dof_index> &renumbering, 1424 std::vector<types::global_dof_index> &block_data, 1425 bool return_start_indices) 1426 { 1427 Assert(renumbering.size() == element.n_dofs_per_cell(), 1428 ExcDimensionMismatch(renumbering.size(), element.n_dofs_per_cell())); 1429 Assert(block_data.size() == element.n_blocks(), 1430 ExcDimensionMismatch(block_data.size(), element.n_blocks())); 1431 1432 types::global_dof_index k = 0; 1433 unsigned int count = 0; 1434 for (unsigned int b = 0; b < element.n_base_elements(); ++b) 1435 for (unsigned int m = 0; m < element.element_multiplicity(b); ++m) 1436 { 1437 block_data[count++] = (return_start_indices) ? 1438 k : 1439 (element.base_element(b).n_dofs_per_cell()); 1440 k += element.base_element(b).n_dofs_per_cell(); 1441 } 1442 Assert(count == element.n_blocks(), ExcInternalError()); 1443 1444 std::vector<types::global_dof_index> start_indices(block_data.size()); 1445 k = 0; 1446 for (unsigned int i = 0; i < block_data.size(); ++i) 1447 if (return_start_indices) 1448 start_indices[i] = block_data[i]; 1449 else 1450 { 1451 start_indices[i] = k; 1452 k += block_data[i]; 1453 } 1454 1455 for (unsigned int i = 0; i < element.n_dofs_per_cell(); ++i) 1456 { 1457 std::pair<unsigned int, types::global_dof_index> indices = 1458 element.system_to_block_index(i); 1459 renumbering[i] = start_indices[indices.first] + indices.second; 1460 } 1461 } 1462 1463 1464 1465 template <int dim, typename number, int spacedim> 1466 void get_interpolation_matrix(const FiniteElement<dim,spacedim> & fe1,const FiniteElement<dim,spacedim> & fe2,FullMatrix<number> & interpolation_matrix)1467 get_interpolation_matrix(const FiniteElement<dim, spacedim> &fe1, 1468 const FiniteElement<dim, spacedim> &fe2, 1469 FullMatrix<number> &interpolation_matrix) 1470 { 1471 Assert(fe1.n_components() == fe2.n_components(), 1472 ExcDimensionMismatch(fe1.n_components(), fe2.n_components())); 1473 Assert(interpolation_matrix.m() == fe2.n_dofs_per_cell() && 1474 interpolation_matrix.n() == fe1.n_dofs_per_cell(), 1475 ExcMatrixDimensionMismatch(interpolation_matrix.m(), 1476 interpolation_matrix.n(), 1477 fe2.n_dofs_per_cell(), 1478 fe1.n_dofs_per_cell())); 1479 1480 // first try the easy way: maybe the FE wants to implement things itself: 1481 try 1482 { 1483 internal::FEToolsGetInterpolationMatrixHelper::gim_forwarder( 1484 fe1, fe2, interpolation_matrix); 1485 return; 1486 } 1487 catch ( 1488 typename FiniteElement<dim, spacedim>::ExcInterpolationNotImplemented &) 1489 { 1490 // too bad.... 1491 } 1492 1493 // uh, so this was not the 1494 // case. hm. then do it the hard 1495 // way. note that this will only 1496 // work if the element is 1497 // primitive, so check this first 1498 Assert(fe1.is_primitive() == true, ExcFENotPrimitive()); 1499 Assert(fe2.is_primitive() == true, ExcFENotPrimitive()); 1500 1501 // Initialize FEValues for fe1 at 1502 // the unit support points of the 1503 // fe2 element. 1504 const std::vector<Point<dim>> &fe2_support_points = 1505 fe2.get_unit_support_points(); 1506 1507 Assert(fe2_support_points.size() == fe2.n_dofs_per_cell(), 1508 (typename FiniteElement<dim, spacedim>::ExcFEHasNoSupportPoints())); 1509 1510 for (unsigned int i = 0; i < fe2.n_dofs_per_cell(); ++i) 1511 { 1512 const unsigned int i1 = fe2.system_to_component_index(i).first; 1513 for (unsigned int j = 0; j < fe1.n_dofs_per_cell(); ++j) 1514 { 1515 const unsigned int j1 = fe1.system_to_component_index(j).first; 1516 if (i1 == j1) 1517 interpolation_matrix(i, j) = 1518 fe1.shape_value(j, fe2_support_points[i]); 1519 else 1520 interpolation_matrix(i, j) = 0.; 1521 } 1522 } 1523 } 1524 1525 1526 1527 template <int dim, typename number, int spacedim> 1528 void get_back_interpolation_matrix(const FiniteElement<dim,spacedim> & fe1,const FiniteElement<dim,spacedim> & fe2,FullMatrix<number> & interpolation_matrix)1529 get_back_interpolation_matrix(const FiniteElement<dim, spacedim> &fe1, 1530 const FiniteElement<dim, spacedim> &fe2, 1531 FullMatrix<number> &interpolation_matrix) 1532 { 1533 Assert(fe1.n_components() == fe2.n_components(), 1534 ExcDimensionMismatch(fe1.n_components(), fe2.n_components())); 1535 Assert(interpolation_matrix.m() == fe1.n_dofs_per_cell() && 1536 interpolation_matrix.n() == fe1.n_dofs_per_cell(), 1537 ExcMatrixDimensionMismatch(interpolation_matrix.m(), 1538 interpolation_matrix.n(), 1539 fe1.n_dofs_per_cell(), 1540 fe1.n_dofs_per_cell())); 1541 1542 FullMatrix<number> first_matrix(fe2.n_dofs_per_cell(), 1543 fe1.n_dofs_per_cell()); 1544 FullMatrix<number> second_matrix(fe1.n_dofs_per_cell(), 1545 fe2.n_dofs_per_cell()); 1546 1547 get_interpolation_matrix(fe1, fe2, first_matrix); 1548 get_interpolation_matrix(fe2, fe1, second_matrix); 1549 1550 // int_matrix=second_matrix*first_matrix 1551 second_matrix.mmult(interpolation_matrix, first_matrix); 1552 } 1553 1554 1555 1556 template <int dim, typename number, int spacedim> 1557 void get_interpolation_difference_matrix(const FiniteElement<dim,spacedim> & fe1,const FiniteElement<dim,spacedim> & fe2,FullMatrix<number> & difference_matrix)1558 get_interpolation_difference_matrix(const FiniteElement<dim, spacedim> &fe1, 1559 const FiniteElement<dim, spacedim> &fe2, 1560 FullMatrix<number> &difference_matrix) 1561 { 1562 Assert(fe1.n_components() == fe2.n_components(), 1563 ExcDimensionMismatch(fe1.n_components(), fe2.n_components())); 1564 Assert(difference_matrix.m() == fe1.n_dofs_per_cell() && 1565 difference_matrix.n() == fe1.n_dofs_per_cell(), 1566 ExcMatrixDimensionMismatch(difference_matrix.m(), 1567 difference_matrix.n(), 1568 fe1.n_dofs_per_cell(), 1569 fe1.n_dofs_per_cell())); 1570 1571 FullMatrix<number> interpolation_matrix(fe1.n_dofs_per_cell()); 1572 get_back_interpolation_matrix(fe1, fe2, interpolation_matrix); 1573 1574 // compute difference 1575 difference_matrix = IdentityMatrix(fe1.n_dofs_per_cell()); 1576 difference_matrix.add(-1, interpolation_matrix); 1577 } 1578 1579 1580 1581 template <int dim, typename number, int spacedim> 1582 void get_projection_matrix(const FiniteElement<dim,spacedim> & fe1,const FiniteElement<dim,spacedim> & fe2,FullMatrix<number> & matrix)1583 get_projection_matrix(const FiniteElement<dim, spacedim> &fe1, 1584 const FiniteElement<dim, spacedim> &fe2, 1585 FullMatrix<number> & matrix) 1586 { 1587 Assert(fe1.n_components() == 1, ExcNotImplemented()); 1588 Assert(fe1.n_components() == fe2.n_components(), 1589 ExcDimensionMismatch(fe1.n_components(), fe2.n_components())); 1590 Assert(matrix.m() == fe2.n_dofs_per_cell() && 1591 matrix.n() == fe1.n_dofs_per_cell(), 1592 ExcMatrixDimensionMismatch(matrix.m(), 1593 matrix.n(), 1594 fe2.n_dofs_per_cell(), 1595 fe1.n_dofs_per_cell())); 1596 matrix = 0; 1597 1598 const unsigned int n1 = fe1.n_dofs_per_cell(); 1599 const unsigned int n2 = fe2.n_dofs_per_cell(); 1600 1601 // First, create a local mass matrix for the unit cell 1602 Triangulation<dim, spacedim> tr; 1603 GridGenerator::hyper_cube(tr); 1604 1605 // Choose a Gauss quadrature rule that is exact up to degree 2n-1 1606 const unsigned int degree = 1607 std::max(fe1.tensor_degree(), fe2.tensor_degree()); 1608 Assert(degree != numbers::invalid_unsigned_int, ExcNotImplemented()); 1609 const QGauss<dim> quadrature(degree + 1); 1610 1611 // Set up FEValues. 1612 const UpdateFlags flags = 1613 update_values | update_quadrature_points | update_JxW_values; 1614 FEValues<dim> val1(fe1, quadrature, update_values); 1615 val1.reinit(tr.begin_active()); 1616 1617 FEValues<dim> val2(fe2, quadrature, flags); 1618 val2.reinit(tr.begin_active()); 1619 1620 // Integrate and invert mass matrix. This happens in the target space 1621 FullMatrix<double> mass(n2, n2); 1622 1623 for (unsigned int k = 0; k < quadrature.size(); ++k) 1624 { 1625 const double dx = val2.JxW(k); 1626 for (unsigned int i = 0; i < n2; ++i) 1627 { 1628 const double v = val2.shape_value(i, k); 1629 for (unsigned int j = 0; j < n2; ++j) 1630 mass(i, j) += v * val2.shape_value(j, k) * dx; 1631 } 1632 } 1633 // Invert the matrix. Gauss-Jordan should be sufficient since we expect 1634 // the mass matrix to be well-conditioned 1635 mass.gauss_jordan(); 1636 1637 // Now, test every function of fe1 with test functions of fe2 and 1638 // compute the projection of each unit vector. 1639 Vector<double> b(n2); 1640 Vector<double> x(n2); 1641 1642 for (unsigned int j = 0; j < n1; ++j) 1643 { 1644 b = 0.; 1645 for (unsigned int i = 0; i < n2; ++i) 1646 for (unsigned int k = 0; k < quadrature.size(); ++k) 1647 { 1648 const double dx = val2.JxW(k); 1649 const double u = val1.shape_value(j, k); 1650 const double v = val2.shape_value(i, k); 1651 b(i) += u * v * dx; 1652 } 1653 1654 // Multiply by the inverse 1655 mass.vmult(x, b); 1656 for (unsigned int i = 0; i < n2; ++i) 1657 matrix(i, j) = x(i); 1658 } 1659 } 1660 1661 1662 1663 template <int dim, int spacedim> 1664 FullMatrix<double> compute_node_matrix(const FiniteElement<dim,spacedim> & fe)1665 compute_node_matrix(const FiniteElement<dim, spacedim> &fe) 1666 { 1667 const unsigned int n_dofs = fe.n_dofs_per_cell(); 1668 1669 FullMatrix<double> N(n_dofs, n_dofs); 1670 1671 Assert(fe.has_generalized_support_points(), ExcNotInitialized()); 1672 Assert(fe.n_components() == dim, ExcNotImplemented()); 1673 1674 const std::vector<Point<dim>> &points = fe.get_generalized_support_points(); 1675 1676 // We need the values of the polynomials in all generalized support 1677 // points. This function specifically works for the case where shape 1678 // functions have 'dim' vector components, so allocate that much space 1679 std::vector<Vector<double>> support_point_values(points.size(), 1680 Vector<double>(dim)); 1681 1682 // In this vector, we store the 1683 // result of the interpolation 1684 std::vector<double> nodal_values(n_dofs); 1685 1686 // Get the values of each shape function in turn. Remember that these 1687 // are the 'raw' shape functions (i.e., where the element has not yet 1688 // computed the expansion coefficients with regard to the basis 1689 // provided by the polynomial space). 1690 for (unsigned int i = 0; i < n_dofs; ++i) 1691 { 1692 // get the values of the current set of shape functions 1693 // at the generalized support points 1694 for (unsigned int k = 0; k < points.size(); ++k) 1695 for (unsigned int d = 0; d < dim; ++d) 1696 { 1697 support_point_values[k][d] = 1698 fe.shape_value_component(i, points[k], d); 1699 Assert(numbers::is_finite(support_point_values[k][d]), 1700 ExcInternalError()); 1701 } 1702 1703 fe.convert_generalized_support_point_values_to_dof_values( 1704 support_point_values, nodal_values); 1705 1706 // Enter the interpolated dofs into the matrix 1707 for (unsigned int j = 0; j < n_dofs; ++j) 1708 { 1709 N(j, i) = nodal_values[j]; 1710 Assert(numbers::is_finite(nodal_values[j]), ExcInternalError()); 1711 } 1712 } 1713 1714 return N; 1715 } 1716 1717 1718 1719 namespace internal 1720 { 1721 namespace FEToolsComputeEmbeddingMatricesHelper 1722 { 1723 template <int dim, typename number, int spacedim> 1724 void compute_embedding_for_shape_function(const unsigned int i,const FiniteElement<dim,spacedim> & fe,const FEValues<dim,spacedim> & coarse,const Householder<double> & H,FullMatrix<number> & this_matrix,const double threshold)1725 compute_embedding_for_shape_function( 1726 const unsigned int i, 1727 const FiniteElement<dim, spacedim> &fe, 1728 const FEValues<dim, spacedim> & coarse, 1729 const Householder<double> & H, 1730 FullMatrix<number> & this_matrix, 1731 const double threshold) 1732 { 1733 const unsigned int n = fe.n_dofs_per_cell(); 1734 const unsigned int nd = fe.n_components(); 1735 const unsigned int nq = coarse.n_quadrature_points; 1736 1737 Vector<number> v_coarse(nq * nd); 1738 Vector<number> v_fine(n); 1739 1740 // The right hand side of 1741 // the least squares 1742 // problem consists of the 1743 // function values of the 1744 // coarse grid function in 1745 // each quadrature point. 1746 if (fe.is_primitive()) 1747 { 1748 const unsigned int d = fe.system_to_component_index(i).first; 1749 const double * phi_i = &coarse.shape_value(i, 0); 1750 1751 for (unsigned int k = 0; k < nq; ++k) 1752 v_coarse(k * nd + d) = phi_i[k]; 1753 } 1754 1755 else 1756 for (unsigned int d = 0; d < nd; ++d) 1757 for (unsigned int k = 0; k < nq; ++k) 1758 v_coarse(k * nd + d) = coarse.shape_value_component(i, k, d); 1759 1760 // solve the least squares 1761 // problem. 1762 const double result = H.least_squares(v_fine, v_coarse); 1763 Assert(result <= threshold, FETools::ExcLeastSquaresError(result)); 1764 // Avoid warnings in release mode 1765 (void)result; 1766 (void)threshold; 1767 1768 // Copy into the result 1769 // matrix. Since the matrix 1770 // maps a coarse grid 1771 // function to a fine grid 1772 // function, the columns 1773 // are fine grid. 1774 for (unsigned int j = 0; j < n; ++j) 1775 this_matrix(j, i) = v_fine(j); 1776 } 1777 1778 1779 1780 template <int dim, typename number, int spacedim> 1781 void compute_embedding_matrices_for_refinement_case(const FiniteElement<dim,spacedim> & fe,std::vector<FullMatrix<number>> & matrices,const unsigned int ref_case,const double threshold)1782 compute_embedding_matrices_for_refinement_case( 1783 const FiniteElement<dim, spacedim> &fe, 1784 std::vector<FullMatrix<number>> & matrices, 1785 const unsigned int ref_case, 1786 const double threshold) 1787 { 1788 const unsigned int n = fe.n_dofs_per_cell(); 1789 const unsigned int nc = 1790 GeometryInfo<dim>::n_children(RefinementCase<dim>(ref_case)); 1791 for (unsigned int i = 0; i < nc; ++i) 1792 { 1793 Assert(matrices[i].n() == n, 1794 ExcDimensionMismatch(matrices[i].n(), n)); 1795 Assert(matrices[i].m() == n, 1796 ExcDimensionMismatch(matrices[i].m(), n)); 1797 } 1798 1799 // Set up meshes, one with a single 1800 // reference cell and refine it once 1801 Triangulation<dim, spacedim> tria; 1802 GridGenerator::hyper_cube(tria, 0, 1); 1803 tria.begin_active()->set_refine_flag(RefinementCase<dim>(ref_case)); 1804 tria.execute_coarsening_and_refinement(); 1805 1806 const unsigned int degree = fe.degree; 1807 QGauss<dim> q_fine(degree + 1); 1808 const unsigned int nq = q_fine.size(); 1809 1810 FEValues<dim, spacedim> fine(fe, 1811 q_fine, 1812 update_quadrature_points | 1813 update_JxW_values | update_values); 1814 1815 // We search for the polynomial on 1816 // the small cell, being equal to 1817 // the coarse polynomial in all 1818 // quadrature points. 1819 1820 // First build the matrix for this 1821 // least squares problem. This 1822 // contains the values of the fine 1823 // cell polynomials in the fine 1824 // cell grid points. 1825 1826 // This matrix is the same for all 1827 // children. 1828 fine.reinit(tria.begin_active()); 1829 const unsigned int nd = fe.n_components(); 1830 FullMatrix<number> A(nq * nd, n); 1831 1832 for (unsigned int j = 0; j < n; ++j) 1833 for (unsigned int d = 0; d < nd; ++d) 1834 for (unsigned int k = 0; k < nq; ++k) 1835 A(k * nd + d, j) = fine.shape_value_component(j, k, d); 1836 1837 Householder<double> H(A); 1838 1839 Threads::TaskGroup<void> task_group; 1840 1841 for (const auto &fine_cell : tria.active_cell_iterators()) 1842 { 1843 fine.reinit(fine_cell); 1844 1845 // evaluate on the coarse cell (which 1846 // is the first -- inactive -- cell on 1847 // the lowest level of the 1848 // triangulation we have created) 1849 const std::vector<Point<spacedim>> &q_points_fine = 1850 fine.get_quadrature_points(); 1851 std::vector<Point<dim>> q_points_coarse(q_points_fine.size()); 1852 for (unsigned int i = 0; i < q_points_fine.size(); ++i) 1853 for (unsigned int j = 0; j < dim; ++j) 1854 q_points_coarse[i](j) = q_points_fine[i](j); 1855 const Quadrature<dim> q_coarse(q_points_coarse, 1856 fine.get_JxW_values()); 1857 FEValues<dim, spacedim> coarse(fe, q_coarse, update_values); 1858 1859 coarse.reinit(tria.begin(0)); 1860 1861 FullMatrix<double> &this_matrix = 1862 matrices[fine_cell->active_cell_index()]; 1863 1864 // Compute this once for each 1865 // coarse grid basis function. can 1866 // spawn subtasks if n is 1867 // sufficiently large so that there 1868 // are more than about 5000 1869 // operations in the inner loop 1870 // (which is basically const * n^2 1871 // operations). 1872 if (n > 30) 1873 { 1874 for (unsigned int i = 0; i < n; ++i) 1875 { 1876 task_group += Threads::new_task( 1877 &compute_embedding_for_shape_function<dim, 1878 number, 1879 spacedim>, 1880 i, 1881 fe, 1882 coarse, 1883 H, 1884 this_matrix, 1885 threshold); 1886 } 1887 task_group.join_all(); 1888 } 1889 else 1890 { 1891 for (unsigned int i = 0; i < n; ++i) 1892 { 1893 compute_embedding_for_shape_function<dim, number, spacedim>( 1894 i, fe, coarse, H, this_matrix, threshold); 1895 } 1896 } 1897 1898 // Remove small entries from 1899 // the matrix 1900 for (unsigned int i = 0; i < this_matrix.m(); ++i) 1901 for (unsigned int j = 0; j < this_matrix.n(); ++j) 1902 if (std::fabs(this_matrix(i, j)) < 1e-12) 1903 this_matrix(i, j) = 0.; 1904 } 1905 } 1906 } // namespace FEToolsComputeEmbeddingMatricesHelper 1907 } // namespace internal 1908 1909 1910 1911 template <int dim, typename number, int spacedim> 1912 void compute_embedding_matrices(const FiniteElement<dim,spacedim> & fe,std::vector<std::vector<FullMatrix<number>>> & matrices,const bool isotropic_only,const double threshold)1913 compute_embedding_matrices(const FiniteElement<dim, spacedim> &fe, 1914 std::vector<std::vector<FullMatrix<number>> 1915 1916 > & matrices, 1917 const bool isotropic_only, 1918 const double threshold) 1919 { 1920 Threads::TaskGroup<void> task_group; 1921 1922 // loop over all possible refinement cases 1923 unsigned int ref_case = (isotropic_only) ? 1924 RefinementCase<dim>::isotropic_refinement : 1925 RefinementCase<dim>::cut_x; 1926 1927 for (; ref_case <= RefinementCase<dim>::isotropic_refinement; ++ref_case) 1928 task_group += Threads::new_task( 1929 &internal::FEToolsComputeEmbeddingMatricesHelper:: 1930 compute_embedding_matrices_for_refinement_case<dim, number, spacedim>, 1931 fe, 1932 matrices[ref_case - 1], 1933 ref_case, 1934 threshold); 1935 1936 task_group.join_all(); 1937 } 1938 1939 1940 template <int dim, typename number, int spacedim> 1941 void compute_face_embedding_matrices(const FiniteElement<dim,spacedim> & fe,FullMatrix<number> (& matrices)[GeometryInfo<dim>::max_children_per_face],const unsigned int face_coarse,const unsigned int face_fine,const double threshold)1942 compute_face_embedding_matrices( 1943 const FiniteElement<dim, spacedim> &fe, 1944 FullMatrix<number> (&matrices)[GeometryInfo<dim>::max_children_per_face], 1945 const unsigned int face_coarse, 1946 const unsigned int face_fine, 1947 const double threshold) 1948 { 1949 const unsigned int face_no = face_coarse; 1950 1951 Assert(face_coarse == 0, ExcNotImplemented()); 1952 Assert(face_fine == 0, ExcNotImplemented()); 1953 1954 const unsigned int nc = GeometryInfo<dim>::max_children_per_face; 1955 const unsigned int n = fe.n_dofs_per_face(face_no); 1956 const unsigned int nd = fe.n_components(); 1957 const unsigned int degree = fe.degree; 1958 1959 const bool normal = fe.conforms(FiniteElementData<dim>::Hdiv); 1960 const bool tangential = fe.conforms(FiniteElementData<dim>::Hcurl); 1961 1962 for (unsigned int i = 0; i < nc; ++i) 1963 { 1964 Assert(matrices[i].n() == n, ExcDimensionMismatch(matrices[i].n(), n)); 1965 Assert(matrices[i].m() == n, ExcDimensionMismatch(matrices[i].m(), n)); 1966 } 1967 1968 // In order to make the loops below 1969 // simpler, we introduce vectors 1970 // containing for indices 0-n the 1971 // number of the corresponding 1972 // shape value on the cell. 1973 std::vector<unsigned int> face_c_dofs(n); 1974 std::vector<unsigned int> face_f_dofs(n); 1975 { 1976 unsigned int face_dof = 0; 1977 for (unsigned int i = 0; 1978 i < ReferenceCell::internal::Info::get_face(fe.reference_cell_type(), 1979 face_no) 1980 .n_vertices(); 1981 ++i) 1982 { 1983 const unsigned int offset_c = 1984 GeometryInfo<dim>::face_to_cell_vertices(face_coarse, i) * 1985 fe.n_dofs_per_vertex(); 1986 const unsigned int offset_f = 1987 GeometryInfo<dim>::face_to_cell_vertices(face_fine, i) * 1988 fe.n_dofs_per_vertex(); 1989 for (unsigned int j = 0; j < fe.n_dofs_per_vertex(); ++j) 1990 { 1991 face_c_dofs[face_dof] = offset_c + j; 1992 face_f_dofs[face_dof] = offset_f + j; 1993 ++face_dof; 1994 } 1995 } 1996 1997 for (unsigned int i = 1; i <= ReferenceCell::internal::Info::get_face( 1998 fe.reference_cell_type(), face_no) 1999 .n_lines(); 2000 ++i) 2001 { 2002 const unsigned int offset_c = 2003 fe.get_first_line_index() + 2004 GeometryInfo<dim>::face_to_cell_lines(face_coarse, i - 1) * 2005 fe.n_dofs_per_line(); 2006 const unsigned int offset_f = 2007 fe.get_first_line_index() + 2008 GeometryInfo<dim>::face_to_cell_lines(face_fine, i - 1) * 2009 fe.n_dofs_per_line(); 2010 for (unsigned int j = 0; j < fe.n_dofs_per_line(); ++j) 2011 { 2012 face_c_dofs[face_dof] = offset_c + j; 2013 face_f_dofs[face_dof] = offset_f + j; 2014 ++face_dof; 2015 } 2016 } 2017 2018 if (dim == 3) 2019 { 2020 const unsigned int offset_c = fe.get_first_quad_index(face_coarse); 2021 const unsigned int offset_f = fe.get_first_quad_index(face_fine); 2022 for (unsigned int j = 0; j < fe.n_dofs_per_quad(face_no); ++j) 2023 { 2024 face_c_dofs[face_dof] = offset_c + j; 2025 face_f_dofs[face_dof] = offset_f + j; 2026 ++face_dof; 2027 } 2028 } 2029 Assert(face_dof == fe.n_dofs_per_face(face_no), ExcInternalError()); 2030 } 2031 2032 // Set up meshes, one with a single 2033 // reference cell and refine it once 2034 Triangulation<dim, spacedim> tria; 2035 GridGenerator::hyper_cube(tria, 0, 1); 2036 tria.refine_global(1); 2037 MappingCartesian<dim> mapping; 2038 2039 // Setup quadrature and FEValues 2040 // for a face. We cannot use 2041 // FEFaceValues and 2042 // FESubfaceValues because of 2043 // some nifty handling of 2044 // refinement cases. Guido stops 2045 // disliking and instead starts 2046 // hating the anisotropic implementation 2047 QGauss<dim - 1> q_gauss(degree + 1); 2048 const Quadrature<dim> q_fine = 2049 QProjector<dim>::project_to_face(fe.reference_cell_type(), 2050 q_gauss, 2051 face_fine); 2052 const unsigned int nq = q_fine.size(); 2053 2054 FEValues<dim> fine(mapping, 2055 fe, 2056 q_fine, 2057 update_quadrature_points | update_JxW_values | 2058 update_values); 2059 2060 // We search for the polynomial on 2061 // the small cell, being equal to 2062 // the coarse polynomial in all 2063 // quadrature points. 2064 2065 // First build the matrix for this 2066 // least squares problem. This 2067 // contains the values of the fine 2068 // cell polynomials in the fine 2069 // cell grid points. 2070 2071 // This matrix is the same for all 2072 // children. 2073 fine.reinit(tria.begin_active()); 2074 FullMatrix<number> A(nq * nd, n); 2075 for (unsigned int j = 0; j < n; ++j) 2076 for (unsigned int k = 0; k < nq; ++k) 2077 if (nd != dim) 2078 for (unsigned int d = 0; d < nd; ++d) 2079 A(k * nd + d, j) = fine.shape_value_component(face_f_dofs[j], k, d); 2080 else 2081 { 2082 if (normal) 2083 A(k * nd, j) = fine.shape_value_component(face_f_dofs[j], k, 0); 2084 if (tangential) 2085 for (unsigned int d = 1; d < dim; ++d) 2086 A(k * nd + d, j) = 2087 fine.shape_value_component(face_f_dofs[j], k, d); 2088 } 2089 2090 Householder<double> H(A); 2091 2092 Vector<number> v_coarse(nq * nd); 2093 Vector<number> v_fine(n); 2094 2095 2096 for (unsigned int cell_number = 0; 2097 cell_number < GeometryInfo<dim>::max_children_per_face; 2098 ++cell_number) 2099 { 2100 const Quadrature<dim> q_coarse = QProjector<dim>::project_to_subface( 2101 fe.reference_cell_type(), q_gauss, face_coarse, cell_number); 2102 FEValues<dim> coarse(mapping, fe, q_coarse, update_values); 2103 2104 typename Triangulation<dim, spacedim>::active_cell_iterator fine_cell = 2105 tria.begin(0)->child(GeometryInfo<dim>::child_cell_on_face( 2106 tria.begin(0)->refinement_case(), face_coarse, cell_number)); 2107 fine.reinit(fine_cell); 2108 coarse.reinit(tria.begin(0)); 2109 2110 FullMatrix<double> &this_matrix = matrices[cell_number]; 2111 2112 // Compute this once for each 2113 // coarse grid basis function 2114 for (unsigned int i = 0; i < n; ++i) 2115 { 2116 // The right hand side of 2117 // the least squares 2118 // problem consists of the 2119 // function values of the 2120 // coarse grid function in 2121 // each quadrature point. 2122 for (unsigned int k = 0; k < nq; ++k) 2123 if (nd != dim) 2124 for (unsigned int d = 0; d < nd; ++d) 2125 v_coarse(k * nd + d) = 2126 coarse.shape_value_component(face_c_dofs[i], k, d); 2127 else 2128 { 2129 if (normal) 2130 v_coarse(k * nd) = 2131 coarse.shape_value_component(face_c_dofs[i], k, 0); 2132 if (tangential) 2133 for (unsigned int d = 1; d < dim; ++d) 2134 v_coarse(k * nd + d) = 2135 coarse.shape_value_component(face_c_dofs[i], k, d); 2136 } 2137 // solve the least squares 2138 // problem. 2139 const double result = H.least_squares(v_fine, v_coarse); 2140 Assert(result <= threshold, FETools::ExcLeastSquaresError(result)); 2141 // Avoid compiler warnings in Release mode 2142 (void)result; 2143 (void)threshold; 2144 2145 // Copy into the result 2146 // matrix. Since the matrix 2147 // maps a coarse grid 2148 // function to a fine grid 2149 // function, the columns 2150 // are fine grid. 2151 for (unsigned int j = 0; j < n; ++j) 2152 this_matrix(j, i) = v_fine(j); 2153 } 2154 // Remove small entries from 2155 // the matrix 2156 for (unsigned int i = 0; i < this_matrix.m(); ++i) 2157 for (unsigned int j = 0; j < this_matrix.n(); ++j) 2158 if (std::fabs(this_matrix(i, j)) < 1e-12) 2159 this_matrix(i, j) = 0.; 2160 } 2161 } 2162 2163 2164 template <int dim, typename number, int spacedim> 2165 void compute_projection_matrices(const FiniteElement<dim,spacedim> & fe,std::vector<std::vector<FullMatrix<number>>> & matrices,const bool isotropic_only)2166 compute_projection_matrices(const FiniteElement<dim, spacedim> &fe, 2167 std::vector<std::vector<FullMatrix<number>> 2168 2169 > & matrices, 2170 const bool isotropic_only) 2171 { 2172 const unsigned int n = fe.n_dofs_per_cell(); 2173 const unsigned int nd = fe.n_components(); 2174 const unsigned int degree = fe.degree; 2175 2176 // prepare FEValues, quadrature etc on 2177 // coarse cell 2178 QGauss<dim> q_fine(degree + 1); 2179 const unsigned int nq = q_fine.size(); 2180 2181 // create mass matrix on coarse cell. 2182 FullMatrix<number> mass(n, n); 2183 { 2184 // set up a triangulation for coarse cell 2185 Triangulation<dim, spacedim> tr; 2186 GridGenerator::hyper_cube(tr, 0, 1); 2187 2188 FEValues<dim, spacedim> coarse(fe, 2189 q_fine, 2190 update_JxW_values | update_values); 2191 2192 typename Triangulation<dim, spacedim>::cell_iterator coarse_cell = 2193 tr.begin(0); 2194 coarse.reinit(coarse_cell); 2195 2196 const std::vector<double> &JxW = coarse.get_JxW_values(); 2197 for (unsigned int i = 0; i < n; ++i) 2198 for (unsigned int j = 0; j < n; ++j) 2199 if (fe. 2200 2201 is_primitive() 2202 2203 ) 2204 { 2205 const double *coarse_i = &coarse.shape_value(i, 0); 2206 const double *coarse_j = &coarse.shape_value(j, 0); 2207 double mass_ij = 0; 2208 for (unsigned int k = 0; k < nq; ++k) 2209 mass_ij += JxW[k] * coarse_i[k] * coarse_j[k]; 2210 mass(i, j) = mass_ij; 2211 } 2212 else 2213 { 2214 double mass_ij = 0; 2215 for (unsigned int d = 0; d < nd; ++d) 2216 for (unsigned int k = 0; k < nq; ++k) 2217 mass_ij += JxW[k] * coarse.shape_value_component(i, k, d) * 2218 coarse.shape_value_component(j, k, d); 2219 mass(i, j) = mass_ij; 2220 } 2221 2222 // invert mass matrix 2223 mass. 2224 2225 gauss_jordan(); 2226 } 2227 2228 2229 auto compute_one_case = 2230 [&fe, &q_fine, n, nd, nq](const unsigned int ref_case, 2231 const FullMatrix<double> &inverse_mass_matrix, 2232 std::vector<FullMatrix<double>> &matrices) { 2233 const unsigned int nc = 2234 GeometryInfo<dim>::n_children(RefinementCase<dim>(ref_case)); 2235 2236 for (unsigned int i = 0; i < nc; ++i) 2237 { 2238 Assert(matrices[i]. 2239 2240 n() 2241 2242 == n, 2243 ExcDimensionMismatch(matrices[i]. 2244 2245 n(), 2246 n 2247 2248 )); 2249 Assert(matrices[i]. 2250 2251 m() 2252 2253 == n, 2254 ExcDimensionMismatch(matrices[i]. 2255 2256 m(), 2257 n 2258 2259 )); 2260 } 2261 2262 // create a respective refinement on the triangulation 2263 Triangulation<dim, spacedim> tr; 2264 GridGenerator::hyper_cube(tr, 0, 1); 2265 tr. 2266 2267 begin_active() 2268 -> 2269 2270 set_refine_flag(RefinementCase<dim>(ref_case)); 2271 tr. 2272 2273 execute_coarsening_and_refinement(); 2274 2275 FEValues<dim, spacedim> fine(StaticMappingQ1<dim, spacedim>::mapping, 2276 fe, 2277 q_fine, 2278 update_quadrature_points | 2279 update_JxW_values | update_values); 2280 2281 typename Triangulation<dim, spacedim>::cell_iterator coarse_cell = 2282 tr.begin(0); 2283 2284 Vector<number> v_coarse(n); 2285 Vector<number> v_fine(n); 2286 2287 for (unsigned int cell_number = 0; cell_number < nc; ++cell_number) 2288 { 2289 FullMatrix<double> &this_matrix = matrices[cell_number]; 2290 2291 // Compute right hand side, which is a fine level basis 2292 // function tested with the coarse level functions. 2293 fine.reinit(coarse_cell->child(cell_number)); 2294 const std::vector<Point<spacedim>> &q_points_fine = 2295 fine.get_quadrature_points(); 2296 std::vector<Point<dim>> q_points_coarse(q_points_fine. 2297 2298 size() 2299 2300 ); 2301 for (unsigned int q = 0; q < q_points_fine. 2302 2303 size(); 2304 2305 ++q) 2306 for (unsigned int j = 0; j < dim; ++j) 2307 q_points_coarse[q](j) = q_points_fine[q](j); 2308 Quadrature<dim> q_coarse(q_points_coarse, fine.get_JxW_values()); 2309 FEValues<dim, spacedim> coarse( 2310 StaticMappingQ1<dim, spacedim>::mapping, 2311 fe, 2312 q_coarse, 2313 update_values); 2314 coarse.reinit(coarse_cell); 2315 2316 // Build RHS 2317 2318 const std::vector<double> &JxW = fine.get_JxW_values(); 2319 2320 // Outer loop over all fine grid shape functions phi_j 2321 for (unsigned int j = 0; j < fe.n_dofs_per_cell(); ++j) 2322 { 2323 for (unsigned int i = 0; i < fe.n_dofs_per_cell(); ++i) 2324 { 2325 if (fe. 2326 2327 is_primitive() 2328 2329 ) 2330 { 2331 const double *coarse_i = &coarse.shape_value(i, 0); 2332 const double *fine_j = &fine.shape_value(j, 0); 2333 2334 double update = 0; 2335 for (unsigned int k = 0; k < nq; ++k) 2336 update += JxW[k] * coarse_i[k] * fine_j[k]; 2337 v_fine(i) = update; 2338 } 2339 else 2340 { 2341 double update = 0; 2342 for (unsigned int d = 0; d < nd; ++d) 2343 for (unsigned int k = 0; k < nq; ++k) 2344 update += JxW[k] * 2345 coarse.shape_value_component(i, k, d) * 2346 fine.shape_value_component(j, k, d); 2347 v_fine(i) = update; 2348 } 2349 } 2350 2351 // RHS ready. Solve system and enter row into matrix 2352 inverse_mass_matrix.vmult(v_coarse, v_fine); 2353 for (unsigned int i = 0; i < fe.n_dofs_per_cell(); ++i) 2354 this_matrix(i, j) = v_coarse(i); 2355 } 2356 2357 // Remove small entries from the matrix 2358 for (unsigned int i = 0; i < this_matrix. 2359 2360 m(); 2361 2362 ++i) 2363 for (unsigned int j = 0; j < this_matrix. 2364 2365 n(); 2366 2367 ++j) 2368 if (std::fabs(this_matrix(i, j)) < 1e-12) 2369 this_matrix(i, j) = 0.; 2370 } 2371 }; 2372 2373 2374 // finally loop over all possible refinement cases 2375 Threads::TaskGroup<> tasks; 2376 unsigned int ref_case = (isotropic_only) ? 2377 RefinementCase<dim>::isotropic_refinement : 2378 RefinementCase<dim>::cut_x; 2379 for (; ref_case <= RefinementCase<dim>::isotropic_refinement; ++ref_case) 2380 tasks += Threads::new_task([&, ref_case]() { 2381 compute_one_case(ref_case, mass, matrices[ref_case - 1]); 2382 }); 2383 2384 tasks. 2385 2386 join_all(); 2387 } 2388 2389 2390 template <int dim, int spacedim> 2391 void add_fe_name(const std::string & parameter_name,const FEFactoryBase<dim,spacedim> * factory)2392 add_fe_name(const std::string & parameter_name, 2393 const FEFactoryBase<dim, spacedim> *factory) 2394 { 2395 // Erase everything after the 2396 // actual class name 2397 std::string name = parameter_name; 2398 unsigned int name_end = name.find_first_not_of(std::string( 2399 "abcdefghijklmnopqrstuvwxyzABCDEFGHIJKLMNOPQRSTUVWXYZ0123456789_")); 2400 if (name_end < name.size()) 2401 name.erase(name_end); 2402 // first make sure that no other 2403 // thread intercepts the 2404 // operation of this function; 2405 // for this, acquire the lock 2406 // until we quit this function 2407 std::lock_guard<std::mutex> lock( 2408 internal::FEToolsAddFENameHelper::fe_name_map_lock); 2409 2410 Assert( 2411 internal::FEToolsAddFENameHelper::get_fe_name_map()[dim][spacedim].find( 2412 name) == 2413 internal::FEToolsAddFENameHelper::get_fe_name_map()[dim][spacedim] 2414 .end(), 2415 ExcMessage("Cannot change existing element in finite element name list")); 2416 2417 // Insert the normalized name into 2418 // the map 2419 internal::FEToolsAddFENameHelper::get_fe_name_map()[dim][spacedim][name] = 2420 std::unique_ptr<const Subscriptor>(factory); 2421 } 2422 2423 2424 2425 namespace internal 2426 { 2427 namespace FEToolsGetFEHelper 2428 { 2429 // TODO: this encapsulates the call to the 2430 // dimension-dependent fe_name_map so that we 2431 // have a unique interface. could be done 2432 // smarter? 2433 template <int dim, int spacedim> 2434 std::unique_ptr<FiniteElement<dim, spacedim>> get_fe_by_name_ext(std::string & name,const std::map<std::string,std::unique_ptr<const Subscriptor>> & fe_name_map)2435 get_fe_by_name_ext( 2436 std::string &name, 2437 const std::map<std::string, std::unique_ptr<const Subscriptor>> 2438 &fe_name_map) 2439 { 2440 // Extract the name of the 2441 // finite element class, which only 2442 // contains characters, numbers and 2443 // underscores. 2444 unsigned int name_end = name.find_first_not_of(std::string( 2445 "abcdefghijklmnopqrstuvwxyzABCDEFGHIJKLMNOPQRSTUVWXYZ0123456789_")); 2446 const std::string name_part(name, 0, name_end); 2447 name.erase(0, name_part.size()); 2448 2449 // now things get a little more 2450 // complicated: FESystem. it's 2451 // more complicated, since we 2452 // have to figure out what the 2453 // base elements are. this can 2454 // only be done recursively 2455 if (name_part == "FESystem") 2456 { 2457 // next we have to get at the 2458 // base elements. start with 2459 // the first. wrap the whole 2460 // block into try-catch to 2461 // make sure we destroy the 2462 // pointers we got from 2463 // recursive calls if one of 2464 // these calls should throw 2465 // an exception 2466 std::vector<std::unique_ptr<const FiniteElement<dim, spacedim>>> 2467 base_fes; 2468 std::vector<unsigned int> base_multiplicities; 2469 2470 // Now, just the [...] 2471 // part should be left. 2472 if (name.size() == 0 || name[0] != '[') 2473 throw(std::string("Invalid first character in ") + name); 2474 do 2475 { 2476 // Erase the 2477 // leading '[' or '-' 2478 name.erase(0, 1); 2479 // Now, the name of the 2480 // first base element is 2481 // first... Let's get it 2482 base_fes.push_back( 2483 get_fe_by_name_ext<dim, spacedim>(name, fe_name_map)); 2484 // next check whether 2485 // FESystem placed a 2486 // multiplicity after 2487 // the element name 2488 if (name[0] == '^') 2489 { 2490 // yes. Delete the '^' 2491 // and read this 2492 // multiplicity 2493 name.erase(0, 1); 2494 2495 const std::pair<int, unsigned int> tmp = 2496 Utilities::get_integer_at_position(name, 0); 2497 name.erase(0, tmp.second); 2498 // add to length, 2499 // including the '^' 2500 base_multiplicities.push_back(tmp.first); 2501 } 2502 else 2503 // no, so 2504 // multiplicity is 2505 // 1 2506 base_multiplicities.push_back(1); 2507 2508 // so that's it for 2509 // this base 2510 // element. base 2511 // elements are 2512 // separated by '-', 2513 // and the list is 2514 // terminated by ']', 2515 // so loop while the 2516 // next character is 2517 // '-' 2518 } 2519 while (name[0] == '-'); 2520 2521 // so we got to the end 2522 // of the '-' separated 2523 // list. make sure that 2524 // we actually had a ']' 2525 // there 2526 if (name.size() == 0 || name[0] != ']') 2527 throw(std::string("Invalid first character in ") + name); 2528 name.erase(0, 1); 2529 // just one more sanity check 2530 Assert((base_fes.size() == base_multiplicities.size()) && 2531 (base_fes.size() > 0), 2532 ExcInternalError()); 2533 2534 // this is a workaround since the constructor for FESystem 2535 // only accepts raw pointers. 2536 std::vector<const FiniteElement<dim, spacedim> *> raw_base_fes; 2537 raw_base_fes.reserve(base_fes.size()); 2538 for (const auto &base_fe : base_fes) 2539 raw_base_fes.push_back(base_fe.get()); 2540 2541 // ok, apparently everything went ok. so generate the composed 2542 // element. and return it. 2543 return std::make_unique<FESystem<dim, spacedim>>( 2544 raw_base_fes, base_multiplicities); 2545 } 2546 else if (name_part == "FE_Nothing") 2547 { 2548 // remove the () from FE_Nothing() 2549 name.erase(0, 2); 2550 2551 // this is a bit of a hack, as 2552 // FE_Nothing does not take a 2553 // degree, but it does take an 2554 // argument, which defaults to 1, 2555 // so this properly returns 2556 // FE_Nothing() 2557 const Subscriptor *ptr = fe_name_map.find(name_part)->second.get(); 2558 const FETools::FEFactoryBase<dim, spacedim> *fef = 2559 dynamic_cast<const FETools::FEFactoryBase<dim, spacedim> *>(ptr); 2560 return fef->get(1); 2561 } 2562 else 2563 { 2564 // Make sure no other thread 2565 // is just adding an element 2566 std::lock_guard<std::mutex> lock( 2567 internal::FEToolsAddFENameHelper::fe_name_map_lock); 2568 AssertThrow(fe_name_map.find(name_part) != fe_name_map.end(), 2569 FETools::ExcInvalidFEName(name)); 2570 2571 // Now, just the (degree) 2572 // or (Quadrature<1>(degree+1)) 2573 // part should be left. 2574 if (name.size() == 0 || name[0] != '(') 2575 throw(std::string("Invalid first character in ") + name); 2576 name.erase(0, 1); 2577 if (name[0] != 'Q') 2578 { 2579 const std::pair<int, unsigned int> tmp = 2580 Utilities::get_integer_at_position(name, 0); 2581 name.erase(0, tmp.second + 1); 2582 const Subscriptor *ptr = 2583 fe_name_map.find(name_part)->second.get(); 2584 const FETools::FEFactoryBase<dim, spacedim> *fef = 2585 dynamic_cast<const FETools::FEFactoryBase<dim, spacedim> *>( 2586 ptr); 2587 return fef->get(tmp.first); 2588 } 2589 else 2590 { 2591 unsigned int position = name.find('('); 2592 const std::string quadrature_name(name, 0, position); 2593 name.erase(0, position + 1); 2594 if (quadrature_name.compare("QGaussLobatto") == 0) 2595 { 2596 const std::pair<int, unsigned int> tmp = 2597 Utilities::get_integer_at_position(name, 0); 2598 // delete "))" 2599 name.erase(0, tmp.second + 2); 2600 const Subscriptor *ptr = 2601 fe_name_map.find(name_part)->second.get(); 2602 const FETools::FEFactoryBase<dim, spacedim> *fef = 2603 dynamic_cast< 2604 const FETools::FEFactoryBase<dim, spacedim> *>(ptr); 2605 return fef->get(QGaussLobatto<1>(tmp.first)); 2606 } 2607 else if (quadrature_name.compare("QGauss") == 0) 2608 { 2609 const std::pair<int, unsigned int> tmp = 2610 Utilities::get_integer_at_position(name, 0); 2611 // delete "))" 2612 name.erase(0, tmp.second + 2); 2613 const Subscriptor *ptr = 2614 fe_name_map.find(name_part)->second.get(); 2615 const FETools::FEFactoryBase<dim, spacedim> *fef = 2616 dynamic_cast< 2617 const FETools::FEFactoryBase<dim, spacedim> *>(ptr); 2618 return fef->get(QGauss<1>(tmp.first)); 2619 } 2620 else if (quadrature_name.compare("QIterated") == 0) 2621 { 2622 // find sub-quadrature 2623 position = name.find('('); 2624 const std::string subquadrature_name(name, 0, position); 2625 AssertThrow(subquadrature_name == "QTrapez" || 2626 subquadrature_name == "QTrapezoid", 2627 ExcNotImplemented( 2628 "Could not detect quadrature of name " + 2629 subquadrature_name)); 2630 // delete "QTrapezoid()," 2631 name.erase(0, position + 3); 2632 const std::pair<int, unsigned int> tmp = 2633 Utilities::get_integer_at_position(name, 0); 2634 // delete "))" 2635 name.erase(0, tmp.second + 2); 2636 const Subscriptor *ptr = 2637 fe_name_map.find(name_part)->second.get(); 2638 const FETools::FEFactoryBase<dim, spacedim> *fef = 2639 dynamic_cast< 2640 const FETools::FEFactoryBase<dim, spacedim> *>(ptr); 2641 return fef->get(QIterated<1>(QTrapezoid<1>(), tmp.first)); 2642 } 2643 else 2644 { 2645 AssertThrow(false, ExcNotImplemented()); 2646 } 2647 } 2648 } 2649 2650 2651 // hm, if we have come thus far, we 2652 // didn't know what to do with the 2653 // string we got. so do as the docs 2654 // say: raise an exception 2655 AssertThrow(false, FETools::ExcInvalidFEName(name)); 2656 2657 // make some compilers happy that 2658 // do not realize that we can't get 2659 // here after throwing 2660 return nullptr; 2661 } 2662 2663 2664 2665 template <int dim, int spacedim> 2666 std::unique_ptr<FiniteElement<dim, spacedim>> get_fe_by_name(std::string & name)2667 get_fe_by_name(std::string &name) 2668 { 2669 return get_fe_by_name_ext<dim, spacedim>( 2670 name, FEToolsAddFENameHelper::get_fe_name_map()[dim][spacedim]); 2671 } 2672 } // namespace FEToolsGetFEHelper 2673 } // namespace internal 2674 2675 2676 2677 template <int dim, int spacedim> 2678 std::unique_ptr<FiniteElement<dim, spacedim>> get_fe_by_name(const std::string & parameter_name)2679 get_fe_by_name(const std::string ¶meter_name) 2680 { 2681 std::string name = Utilities::trim(parameter_name); 2682 std::size_t index = 1; 2683 // remove spaces that are not between two word (things that match the 2684 // regular expression [A-Za-z0-9_]) characters. 2685 while (2 < name.size() && index < name.size() - 1) 2686 { 2687 if (name[index] == ' ' && 2688 (!(std::isalnum(name[index - 1]) || name[index - 1] == '_') || 2689 !(std::isalnum(name[index + 1]) || name[index + 1] == '_'))) 2690 { 2691 name.erase(index, 1); 2692 } 2693 else 2694 { 2695 ++index; 2696 } 2697 } 2698 2699 // Create a version of the name 2700 // string where all template 2701 // parameters are eliminated. 2702 for (unsigned int pos1 = name.find('<'); pos1 < name.size(); 2703 pos1 = name.find('<')) 2704 { 2705 const unsigned int pos2 = name.find('>'); 2706 // If there is only a single 2707 // character between those two, 2708 // it should be 'd' or the number 2709 // representing the dimension. 2710 if (pos2 - pos1 == 2) 2711 { 2712 const char dimchar = '0' + dim; 2713 (void)dimchar; 2714 if (name.at(pos1 + 1) != 'd') 2715 Assert(name.at(pos1 + 1) == dimchar, 2716 ExcInvalidFEDimension(name.at(pos1 + 1), dim)); 2717 } 2718 else 2719 Assert(pos2 - pos1 == 4, ExcInvalidFEName(name)); 2720 2721 // If pos1==pos2, then we are 2722 // probably at the end of the 2723 // string 2724 if (pos2 != pos1) 2725 name.erase(pos1, pos2 - pos1 + 1); 2726 } 2727 // Replace all occurrences of "^dim" 2728 // by "^d" to be handled by the 2729 // next loop 2730 for (unsigned int pos = name.find("^dim"); pos < name.size(); 2731 pos = name.find("^dim")) 2732 name.erase(pos + 2, 2); 2733 2734 // Replace all occurrences of "^d" 2735 // by using the actual dimension 2736 for (unsigned int pos = name.find("^d"); pos < name.size(); 2737 pos = name.find("^d")) 2738 name.at(pos + 1) = '0' + dim; 2739 2740 try 2741 { 2742 auto fe = 2743 internal::FEToolsGetFEHelper::get_fe_by_name<dim, spacedim>(name); 2744 2745 // Make sure the auxiliary function 2746 // ate up all characters of the name. 2747 AssertThrow(name.size() == 0, 2748 ExcInvalidFEName(parameter_name + 2749 std::string(" extra characters after " 2750 "end of name"))); 2751 return fe; 2752 } 2753 catch (const std::string &errline) 2754 { 2755 AssertThrow(false, 2756 ExcInvalidFEName(parameter_name + std::string(" at ") + 2757 errline)); 2758 return nullptr; 2759 } 2760 } 2761 2762 2763 2764 template <int dim, int spacedim> 2765 void compute_projection_from_quadrature_points_matrix(const FiniteElement<dim,spacedim> & fe,const Quadrature<dim> & lhs_quadrature,const Quadrature<dim> & rhs_quadrature,FullMatrix<double> & X)2766 compute_projection_from_quadrature_points_matrix( 2767 const FiniteElement<dim, spacedim> &fe, 2768 const Quadrature<dim> & lhs_quadrature, 2769 const Quadrature<dim> & rhs_quadrature, 2770 FullMatrix<double> & X) 2771 { 2772 Assert(fe.n_components() == 1, ExcNotImplemented()); 2773 2774 // first build the matrices M and Q 2775 // described in the documentation 2776 FullMatrix<double> M(fe.n_dofs_per_cell(), fe.n_dofs_per_cell()); 2777 FullMatrix<double> Q(fe.n_dofs_per_cell(), rhs_quadrature.size()); 2778 2779 for (unsigned int i = 0; i < fe.n_dofs_per_cell(); ++i) 2780 for (unsigned int j = 0; j < fe.n_dofs_per_cell(); ++j) 2781 for (unsigned int q = 0; q < lhs_quadrature.size(); ++q) 2782 M(i, j) += fe.shape_value(i, lhs_quadrature.point(q)) * 2783 fe.shape_value(j, lhs_quadrature.point(q)) * 2784 lhs_quadrature.weight(q); 2785 2786 for (unsigned int i = 0; i < fe.n_dofs_per_cell(); ++i) 2787 for (unsigned int q = 0; q < rhs_quadrature.size(); ++q) 2788 Q(i, q) += 2789 fe.shape_value(i, rhs_quadrature.point(q)) * rhs_quadrature.weight(q); 2790 2791 // then invert M 2792 FullMatrix<double> M_inverse(fe.n_dofs_per_cell(), fe.n_dofs_per_cell()); 2793 M_inverse.invert(M); 2794 2795 // finally compute the result 2796 X.reinit(fe.n_dofs_per_cell(), rhs_quadrature.size()); 2797 M_inverse.mmult(X, Q); 2798 2799 Assert(X.m() == fe.n_dofs_per_cell(), ExcInternalError()); 2800 Assert(X.n() == rhs_quadrature.size(), ExcInternalError()); 2801 } 2802 2803 2804 2805 template <int dim, int spacedim> 2806 void compute_interpolation_to_quadrature_points_matrix(const FiniteElement<dim,spacedim> & fe,const Quadrature<dim> & quadrature,FullMatrix<double> & I_q)2807 compute_interpolation_to_quadrature_points_matrix( 2808 const FiniteElement<dim, spacedim> &fe, 2809 const Quadrature<dim> & quadrature, 2810 FullMatrix<double> & I_q) 2811 { 2812 Assert(fe.n_components() == 1, ExcNotImplemented()); 2813 Assert(I_q.m() == quadrature.size(), ExcMessage("Wrong matrix size")); 2814 Assert(I_q.n() == fe.n_dofs_per_cell(), ExcMessage("Wrong matrix size")); 2815 2816 for (unsigned int q = 0; q < quadrature.size(); ++q) 2817 for (unsigned int i = 0; i < fe.n_dofs_per_cell(); ++i) 2818 I_q(q, i) = fe.shape_value(i, quadrature.point(q)); 2819 } 2820 2821 2822 2823 template <int dim> 2824 void compute_projection_from_quadrature_points(const FullMatrix<double> & projection_matrix,const std::vector<Tensor<1,dim>> & vector_of_tensors_at_qp,std::vector<Tensor<1,dim>> & vector_of_tensors_at_nodes)2825 compute_projection_from_quadrature_points( 2826 const FullMatrix<double> & projection_matrix, 2827 const std::vector<Tensor<1, dim>> &vector_of_tensors_at_qp, 2828 std::vector<Tensor<1, dim>> & vector_of_tensors_at_nodes) 2829 { 2830 // check that the number columns of the projection_matrix 2831 // matches the size of the vector_of_tensors_at_qp 2832 Assert(projection_matrix.n_cols() == vector_of_tensors_at_qp.size(), 2833 ExcDimensionMismatch(projection_matrix.n_cols(), 2834 vector_of_tensors_at_qp.size())); 2835 2836 // check that the number rows of the projection_matrix 2837 // matches the size of the vector_of_tensors_at_nodes 2838 Assert(projection_matrix.n_rows() == vector_of_tensors_at_nodes.size(), 2839 ExcDimensionMismatch(projection_matrix.n_rows(), 2840 vector_of_tensors_at_nodes.size())); 2841 2842 // number of support points (nodes) to project to 2843 const unsigned int n_support_points = projection_matrix.n_rows(); 2844 // number of quadrature points to project from 2845 const unsigned int n_quad_points = projection_matrix.n_cols(); 2846 2847 // component projected to the nodes 2848 Vector<double> component_at_node(n_support_points); 2849 // component at the quadrature point 2850 Vector<double> component_at_qp(n_quad_points); 2851 2852 for (unsigned int ii = 0; ii < dim; ++ii) 2853 { 2854 component_at_qp = 0; 2855 2856 // populate the vector of components at the qps 2857 // from vector_of_tensors_at_qp 2858 // vector_of_tensors_at_qp data is in form: 2859 // columns: 0, 1, ..., dim 2860 // rows: 0,1,...., n_quad_points 2861 // so extract the ii'th column of vector_of_tensors_at_qp 2862 for (unsigned int q = 0; q < n_quad_points; ++q) 2863 { 2864 component_at_qp(q) = vector_of_tensors_at_qp[q][ii]; 2865 } 2866 2867 // project from the qps -> nodes 2868 // component_at_node = projection_matrix_u * component_at_qp 2869 projection_matrix.vmult(component_at_node, component_at_qp); 2870 2871 // rewrite the projection of the components 2872 // back into the vector of tensors 2873 for (unsigned int nn = 0; nn < n_support_points; ++nn) 2874 { 2875 vector_of_tensors_at_nodes[nn][ii] = component_at_node(nn); 2876 } 2877 } 2878 } 2879 2880 2881 2882 template <int dim> 2883 void compute_projection_from_quadrature_points(const FullMatrix<double> & projection_matrix,const std::vector<SymmetricTensor<2,dim>> & vector_of_tensors_at_qp,std::vector<SymmetricTensor<2,dim>> & vector_of_tensors_at_nodes)2884 compute_projection_from_quadrature_points( 2885 const FullMatrix<double> & projection_matrix, 2886 const std::vector<SymmetricTensor<2, dim>> &vector_of_tensors_at_qp, 2887 std::vector<SymmetricTensor<2, dim>> & vector_of_tensors_at_nodes) 2888 { 2889 // check that the number columns of the projection_matrix 2890 // matches the size of the vector_of_tensors_at_qp 2891 Assert(projection_matrix.n_cols() == vector_of_tensors_at_qp.size(), 2892 ExcDimensionMismatch(projection_matrix.n_cols(), 2893 vector_of_tensors_at_qp.size())); 2894 2895 // check that the number rows of the projection_matrix 2896 // matches the size of the vector_of_tensors_at_nodes 2897 Assert(projection_matrix.n_rows() == vector_of_tensors_at_nodes.size(), 2898 ExcDimensionMismatch(projection_matrix.n_rows(), 2899 vector_of_tensors_at_nodes.size())); 2900 2901 // number of support points (nodes) 2902 const unsigned int n_support_points = projection_matrix.n_rows(); 2903 // number of quadrature points to project from 2904 const unsigned int n_quad_points = projection_matrix.n_cols(); 2905 2906 // number of unique entries in a symmetric second-order tensor 2907 const unsigned int n_independent_components = 2908 SymmetricTensor<2, dim>::n_independent_components; 2909 2910 // component projected to the nodes 2911 Vector<double> component_at_node(n_support_points); 2912 // component at the quadrature point 2913 Vector<double> component_at_qp(n_quad_points); 2914 2915 // loop over the number of unique dimensions of the tensor 2916 for (unsigned int ii = 0; ii < n_independent_components; ++ii) 2917 { 2918 component_at_qp = 0; 2919 2920 // row-column entry of tensor corresponding the unrolled index 2921 TableIndices<2> row_column_index = 2922 SymmetricTensor<2, dim>::unrolled_to_component_indices(ii); 2923 const unsigned int row = row_column_index[0]; 2924 const unsigned int column = row_column_index[1]; 2925 2926 // populate the vector of components at the qps 2927 // from vector_of_tensors_at_qp 2928 // vector_of_tensors_at_qp is in form: 2929 // columns: 0, 1, ..., n_independent_components 2930 // rows: 0,1,...., n_quad_points 2931 // so extract the ii'th column of vector_of_tensors_at_qp 2932 for (unsigned int q = 0; q < n_quad_points; ++q) 2933 { 2934 component_at_qp(q) = (vector_of_tensors_at_qp[q])[row][column]; 2935 } 2936 2937 // project from the qps -> nodes 2938 // component_at_node = projection_matrix_u * component_at_qp 2939 projection_matrix.vmult(component_at_node, component_at_qp); 2940 2941 // rewrite the projection of the components back into the vector of 2942 // tensors 2943 for (unsigned int nn = 0; nn < n_support_points; ++nn) 2944 { 2945 (vector_of_tensors_at_nodes[nn])[row][column] = 2946 component_at_node(nn); 2947 } 2948 } 2949 } 2950 2951 2952 2953 template <int dim, int spacedim> 2954 void compute_projection_from_face_quadrature_points_matrix(const FiniteElement<dim,spacedim> & fe,const Quadrature<dim-1> & lhs_quadrature,const Quadrature<dim-1> & rhs_quadrature,const typename DoFHandler<dim,spacedim>::active_cell_iterator & cell,const unsigned int face,FullMatrix<double> & X)2955 compute_projection_from_face_quadrature_points_matrix( 2956 const FiniteElement<dim, spacedim> &fe, 2957 const Quadrature<dim - 1> & lhs_quadrature, 2958 const Quadrature<dim - 1> & rhs_quadrature, 2959 const typename DoFHandler<dim, spacedim>::active_cell_iterator &cell, 2960 const unsigned int face, 2961 FullMatrix<double> & X) 2962 { 2963 Assert(fe.n_components() == 1, ExcNotImplemented()); 2964 Assert(lhs_quadrature.size() > fe.degree, 2965 ExcNotGreaterThan(lhs_quadrature.size(), fe.degree)); 2966 2967 2968 2969 // build the matrices M and Q 2970 // described in the documentation 2971 FullMatrix<double> M(fe.n_dofs_per_cell(), fe.n_dofs_per_cell()); 2972 FullMatrix<double> Q(fe.n_dofs_per_cell(), rhs_quadrature.size()); 2973 2974 { 2975 // need an FEFaceValues object to evaluate shape function 2976 // values on the specified face. 2977 FEFaceValues<dim> fe_face_values(fe, lhs_quadrature, update_values); 2978 fe_face_values.reinit(cell, face); // setup shape_value on this face. 2979 2980 for (unsigned int i = 0; i < fe.n_dofs_per_cell(); ++i) 2981 for (unsigned int j = 0; j < fe.n_dofs_per_cell(); ++j) 2982 for (unsigned int q = 0; q < lhs_quadrature.size(); ++q) 2983 M(i, j) += fe_face_values.shape_value(i, q) * 2984 fe_face_values.shape_value(j, q) * 2985 lhs_quadrature.weight(q); 2986 for (unsigned int i = 0; i < fe.n_dofs_per_cell(); ++i) 2987 { 2988 M(i, i) = (M(i, i) == 0 ? 1 : M(i, i)); 2989 } 2990 } 2991 2992 { 2993 FEFaceValues<dim> fe_face_values(fe, rhs_quadrature, update_values); 2994 fe_face_values.reinit(cell, face); // setup shape_value on this face. 2995 2996 for (unsigned int i = 0; i < fe.n_dofs_per_cell(); ++i) 2997 for (unsigned int q = 0; q < rhs_quadrature.size(); ++q) 2998 Q(i, q) += 2999 fe_face_values.shape_value(i, q) * rhs_quadrature.weight(q); 3000 } 3001 // then invert M 3002 FullMatrix<double> M_inverse(fe.n_dofs_per_cell(), fe.n_dofs_per_cell()); 3003 M_inverse.invert(M); 3004 3005 // finally compute the result 3006 X.reinit(fe.n_dofs_per_cell(), rhs_quadrature.size()); 3007 M_inverse.mmult(X, Q); 3008 3009 Assert(X.m() == fe.n_dofs_per_cell(), ExcInternalError()); 3010 Assert(X.n() == rhs_quadrature.size(), ExcInternalError()); 3011 } 3012 3013 3014 3015 namespace internal 3016 { 3017 namespace FEToolsConvertHelper 3018 { 3019 // Helper functions for 3020 // FETools::convert_generalized_support_point_values_to_dof_values 3021 3022 template <int dim, int spacedim, typename number> 3023 static void convert_helper(const FiniteElement<dim,spacedim> & finite_element,const std::vector<Vector<number>> & support_point_values,std::vector<number> & dof_values)3024 convert_helper(const FiniteElement<dim, spacedim> &finite_element, 3025 const std::vector<Vector<number>> & support_point_values, 3026 std::vector<number> & dof_values) 3027 { 3028 static Threads::ThreadLocalStorage<std::vector<Vector<double>>> 3029 double_support_point_values; 3030 static Threads::ThreadLocalStorage<std::vector<double>> 3031 double_dof_values; 3032 3033 double_support_point_values.get().resize(support_point_values.size()); 3034 double_dof_values.get().resize(dof_values.size()); 3035 3036 for (unsigned int i = 0; i < support_point_values.size(); ++i) 3037 { 3038 double_support_point_values.get()[i].reinit( 3039 finite_element.n_components(), false); 3040 std::copy(std::begin(support_point_values[i]), 3041 std::end(support_point_values[i]), 3042 std::begin(double_support_point_values.get()[i])); 3043 } 3044 3045 finite_element.convert_generalized_support_point_values_to_dof_values( 3046 double_support_point_values.get(), double_dof_values.get()); 3047 3048 std::copy(std::begin(double_dof_values.get()), 3049 std::end(double_dof_values.get()), 3050 std::begin(dof_values)); 3051 } 3052 3053 3054 template <int dim, int spacedim, typename number> 3055 static void convert_helper(const FiniteElement<dim,spacedim> & finite_element,const std::vector<Vector<std::complex<number>>> & support_point_values,std::vector<std::complex<number>> & dof_values)3056 convert_helper( 3057 const FiniteElement<dim, spacedim> & finite_element, 3058 const std::vector<Vector<std::complex<number>>> &support_point_values, 3059 std::vector<std::complex<number>> & dof_values) 3060 { 3061 static Threads::ThreadLocalStorage<std::vector<Vector<double>>> 3062 double_support_point_values_real; 3063 static Threads::ThreadLocalStorage<std::vector<double>> 3064 double_dof_values_real; 3065 static Threads::ThreadLocalStorage<std::vector<Vector<double>>> 3066 double_support_point_values_imag; 3067 static Threads::ThreadLocalStorage<std::vector<double>> 3068 double_dof_values_imag; 3069 3070 double_support_point_values_real.get().resize( 3071 support_point_values.size()); 3072 double_dof_values_real.get().resize(dof_values.size()); 3073 double_support_point_values_imag.get().resize( 3074 support_point_values.size()); 3075 double_dof_values_imag.get().resize(dof_values.size()); 3076 3077 for (unsigned int i = 0; i < support_point_values.size(); ++i) 3078 { 3079 double_support_point_values_real.get()[i].reinit( 3080 finite_element.n_components(), false); 3081 double_support_point_values_imag.get()[i].reinit( 3082 finite_element.n_components(), false); 3083 3084 std::transform( 3085 std::begin(support_point_values[i]), 3086 std::end(support_point_values[i]), 3087 std::begin(double_support_point_values_real.get()[i]), 3088 [](std::complex<number> c) -> double { return c.real(); }); 3089 3090 std::transform( 3091 std::begin(support_point_values[i]), 3092 std::end(support_point_values[i]), 3093 std::begin(double_support_point_values_imag.get()[i]), 3094 [](std::complex<number> c) -> double { return c.imag(); }); 3095 } 3096 3097 finite_element.convert_generalized_support_point_values_to_dof_values( 3098 double_support_point_values_real.get(), double_dof_values_real.get()); 3099 finite_element.convert_generalized_support_point_values_to_dof_values( 3100 double_support_point_values_imag.get(), double_dof_values_imag.get()); 3101 3102 std::transform(std::begin(double_dof_values_real.get()), 3103 std::end(double_dof_values_real.get()), 3104 std::begin(double_dof_values_imag.get()), 3105 std::begin(dof_values), 3106 [](number real, number imag) -> std::complex<number> { 3107 return {real, imag}; 3108 }); 3109 } 3110 3111 3112 template <int dim, int spacedim> 3113 static void convert_helper(const FiniteElement<dim,spacedim> & finite_element,const std::vector<Vector<double>> & support_point_values,std::vector<double> & dof_values)3114 convert_helper(const FiniteElement<dim, spacedim> &finite_element, 3115 const std::vector<Vector<double>> & support_point_values, 3116 std::vector<double> & dof_values) 3117 { 3118 finite_element.convert_generalized_support_point_values_to_dof_values( 3119 support_point_values, dof_values); 3120 } 3121 3122 } // namespace FEToolsConvertHelper 3123 } // namespace internal 3124 3125 3126 3127 template <int dim, int spacedim, typename number> 3128 void convert_generalized_support_point_values_to_dof_values(const FiniteElement<dim,spacedim> & finite_element,const std::vector<Vector<number>> & support_point_values,std::vector<number> & dof_values)3129 convert_generalized_support_point_values_to_dof_values( 3130 const FiniteElement<dim, spacedim> &finite_element, 3131 const std::vector<Vector<number>> & support_point_values, 3132 std::vector<number> & dof_values) 3133 { 3134 AssertDimension(support_point_values.size(), 3135 finite_element.get_generalized_support_points().size()); 3136 AssertDimension(dof_values.size(), finite_element.n_dofs_per_cell()); 3137 3138 internal::FEToolsConvertHelper::convert_helper<dim, spacedim>( 3139 finite_element, support_point_values, dof_values); 3140 } 3141 3142 3143 3144 template <int dim> 3145 std::vector<unsigned int> hierarchic_to_lexicographic_numbering(const unsigned int degree)3146 hierarchic_to_lexicographic_numbering(const unsigned int degree) 3147 { 3148 // number of support points in each direction 3149 const unsigned int n = degree + 1; 3150 3151 const unsigned int dofs_per_cell = Utilities::fixed_power<dim>(n); 3152 3153 std::vector<unsigned int> h2l(dofs_per_cell); 3154 3155 // polynomial degree 3156 const unsigned int dofs_per_line = degree - 1; 3157 3158 // the following lines of code are somewhat odd, due to the way the 3159 // hierarchic numbering is organized. if someone would really want to 3160 // understand these lines, you better draw some pictures where you 3161 // indicate the indices and orders of vertices, lines, etc, along with the 3162 // numbers of the degrees of freedom in hierarchical and lexicographical 3163 // order 3164 switch (dim) 3165 { 3166 case 0: 3167 { 3168 h2l[0] = 0; 3169 break; 3170 } 3171 case 1: 3172 { 3173 h2l[0] = 0; 3174 h2l[1] = dofs_per_cell - 1; 3175 for (unsigned int i = 2; i < dofs_per_cell; ++i) 3176 h2l[i] = i - 1; 3177 3178 break; 3179 } 3180 3181 case 2: 3182 { 3183 unsigned int next_index = 0; 3184 // first the four vertices 3185 h2l[next_index++] = 0; 3186 h2l[next_index++] = n - 1; 3187 h2l[next_index++] = n * (n - 1); 3188 h2l[next_index++] = n * n - 1; 3189 3190 // left line 3191 for (unsigned int i = 0; i < dofs_per_line; ++i) 3192 h2l[next_index++] = (1 + i) * n; 3193 3194 // right line 3195 for (unsigned int i = 0; i < dofs_per_line; ++i) 3196 h2l[next_index++] = (2 + i) * n - 1; 3197 3198 // bottom line 3199 for (unsigned int i = 0; i < dofs_per_line; ++i) 3200 h2l[next_index++] = 1 + i; 3201 3202 // top line 3203 for (unsigned int i = 0; i < dofs_per_line; ++i) 3204 h2l[next_index++] = n * (n - 1) + i + 1; 3205 3206 // inside quad 3207 for (unsigned int i = 0; i < dofs_per_line; ++i) 3208 for (unsigned int j = 0; j < dofs_per_line; ++j) 3209 h2l[next_index++] = n * (i + 1) + j + 1; 3210 3211 Assert(next_index == dofs_per_cell, ExcInternalError()); 3212 3213 break; 3214 } 3215 3216 case 3: 3217 { 3218 unsigned int next_index = 0; 3219 // first the eight vertices 3220 h2l[next_index++] = 0; // 0 3221 h2l[next_index++] = (1) * degree; // 1 3222 h2l[next_index++] = (n)*degree; // 2 3223 h2l[next_index++] = (n + 1) * degree; // 3 3224 h2l[next_index++] = (n * n) * degree; // 4 3225 h2l[next_index++] = (n * n + 1) * degree; // 5 3226 h2l[next_index++] = (n * n + n) * degree; // 6 3227 h2l[next_index++] = (n * n + n + 1) * degree; // 7 3228 3229 // line 0 3230 for (unsigned int i = 0; i < dofs_per_line; ++i) 3231 h2l[next_index++] = (i + 1) * n; 3232 // line 1 3233 for (unsigned int i = 0; i < dofs_per_line; ++i) 3234 h2l[next_index++] = n - 1 + (i + 1) * n; 3235 // line 2 3236 for (unsigned int i = 0; i < dofs_per_line; ++i) 3237 h2l[next_index++] = 1 + i; 3238 // line 3 3239 for (unsigned int i = 0; i < dofs_per_line; ++i) 3240 h2l[next_index++] = 1 + i + n * (n - 1); 3241 3242 // line 4 3243 for (unsigned int i = 0; i < dofs_per_line; ++i) 3244 h2l[next_index++] = (n - 1) * n * n + (i + 1) * n; 3245 // line 5 3246 for (unsigned int i = 0; i < dofs_per_line; ++i) 3247 h2l[next_index++] = (n - 1) * (n * n + 1) + (i + 1) * n; 3248 // line 6 3249 for (unsigned int i = 0; i < dofs_per_line; ++i) 3250 h2l[next_index++] = n * n * (n - 1) + i + 1; 3251 // line 7 3252 for (unsigned int i = 0; i < dofs_per_line; ++i) 3253 h2l[next_index++] = n * n * (n - 1) + i + 1 + n * (n - 1); 3254 3255 // line 8 3256 for (unsigned int i = 0; i < dofs_per_line; ++i) 3257 h2l[next_index++] = (i + 1) * n * n; 3258 // line 9 3259 for (unsigned int i = 0; i < dofs_per_line; ++i) 3260 h2l[next_index++] = n - 1 + (i + 1) * n * n; 3261 // line 10 3262 for (unsigned int i = 0; i < dofs_per_line; ++i) 3263 h2l[next_index++] = (i + 1) * n * n + n * (n - 1); 3264 // line 11 3265 for (unsigned int i = 0; i < dofs_per_line; ++i) 3266 h2l[next_index++] = n - 1 + (i + 1) * n * n + n * (n - 1); 3267 3268 3269 // inside quads 3270 // face 0 3271 for (unsigned int i = 0; i < dofs_per_line; ++i) 3272 for (unsigned int j = 0; j < dofs_per_line; ++j) 3273 h2l[next_index++] = (i + 1) * n * n + n * (j + 1); 3274 // face 1 3275 for (unsigned int i = 0; i < dofs_per_line; ++i) 3276 for (unsigned int j = 0; j < dofs_per_line; ++j) 3277 h2l[next_index++] = (i + 1) * n * n + n - 1 + n * (j + 1); 3278 // face 2, note the orientation! 3279 for (unsigned int i = 0; i < dofs_per_line; ++i) 3280 for (unsigned int j = 0; j < dofs_per_line; ++j) 3281 h2l[next_index++] = (j + 1) * n * n + i + 1; 3282 // face 3, note the orientation! 3283 for (unsigned int i = 0; i < dofs_per_line; ++i) 3284 for (unsigned int j = 0; j < dofs_per_line; ++j) 3285 h2l[next_index++] = (j + 1) * n * n + n * (n - 1) + i + 1; 3286 // face 4 3287 for (unsigned int i = 0; i < dofs_per_line; ++i) 3288 for (unsigned int j = 0; j < dofs_per_line; ++j) 3289 h2l[next_index++] = n * (i + 1) + j + 1; 3290 // face 5 3291 for (unsigned int i = 0; i < dofs_per_line; ++i) 3292 for (unsigned int j = 0; j < dofs_per_line; ++j) 3293 h2l[next_index++] = (n - 1) * n * n + n * (i + 1) + j + 1; 3294 3295 // inside hex 3296 for (unsigned int i = 0; i < dofs_per_line; ++i) 3297 for (unsigned int j = 0; j < dofs_per_line; ++j) 3298 for (unsigned int k = 0; k < dofs_per_line; ++k) 3299 h2l[next_index++] = n * n * (i + 1) + n * (j + 1) + k + 1; 3300 3301 Assert(next_index == dofs_per_cell, ExcInternalError()); 3302 3303 break; 3304 } 3305 3306 default: 3307 Assert(false, ExcNotImplemented()); 3308 } 3309 3310 return h2l; 3311 } 3312 3313 3314 3315 template <int dim> 3316 void hierarchic_to_lexicographic_numbering(const unsigned int degree,std::vector<unsigned int> & h2l)3317 hierarchic_to_lexicographic_numbering(const unsigned int degree, 3318 std::vector<unsigned int> &h2l) 3319 { 3320 AssertDimension(h2l.size(), Utilities::fixed_power<dim>(degree + 1)); 3321 h2l = hierarchic_to_lexicographic_numbering<dim>(degree); 3322 } 3323 3324 3325 3326 template <int dim> 3327 void hierarchic_to_lexicographic_numbering(const FiniteElementData<dim> & fe,std::vector<unsigned int> & h2l)3328 hierarchic_to_lexicographic_numbering(const FiniteElementData<dim> &fe, 3329 std::vector<unsigned int> & h2l) 3330 { 3331 Assert(h2l.size() == fe.n_dofs_per_cell(), 3332 ExcDimensionMismatch(h2l.size(), fe.n_dofs_per_cell())); 3333 hierarchic_to_lexicographic_numbering<dim>(fe.n_dofs_per_line() + 1, h2l); 3334 } 3335 3336 3337 3338 template <int dim> 3339 std::vector<unsigned int> hierarchic_to_lexicographic_numbering(const FiniteElementData<dim> & fe)3340 hierarchic_to_lexicographic_numbering(const FiniteElementData<dim> &fe) 3341 { 3342 Assert(fe.n_components() == 1, ExcInvalidFE()); 3343 return hierarchic_to_lexicographic_numbering<dim>(fe.n_dofs_per_line() + 1); 3344 } 3345 3346 3347 3348 template <int dim> 3349 std::vector<unsigned int> lexicographic_to_hierarchic_numbering(const unsigned int degree)3350 lexicographic_to_hierarchic_numbering(const unsigned int degree) 3351 { 3352 return Utilities::invert_permutation( 3353 hierarchic_to_lexicographic_numbering<dim>(degree)); 3354 } 3355 3356 3357 3358 template <int dim> 3359 void lexicographic_to_hierarchic_numbering(const FiniteElementData<dim> & fe,std::vector<unsigned int> & l2h)3360 lexicographic_to_hierarchic_numbering(const FiniteElementData<dim> &fe, 3361 std::vector<unsigned int> & l2h) 3362 { 3363 l2h = lexicographic_to_hierarchic_numbering(fe); 3364 } 3365 3366 3367 3368 template <int dim> 3369 std::vector<unsigned int> lexicographic_to_hierarchic_numbering(const FiniteElementData<dim> & fe)3370 lexicographic_to_hierarchic_numbering(const FiniteElementData<dim> &fe) 3371 { 3372 return Utilities::invert_permutation( 3373 hierarchic_to_lexicographic_numbering(fe)); 3374 } 3375 3376 } // namespace FETools 3377 3378 3379 DEAL_II_NAMESPACE_CLOSE 3380 3381 #endif /* dealii_fe_tools_templates_H */ 3382