Fix DoFTools::extract_boundary_dof() for simplex meshes

diff --git a/source/dofs/dof_tools.cc b/source/dofs/dof_tools.cc
index 55d820ce8ca409fd18e9360228c2256e5c36e022..30f8bf8980d82ef647cbdc57b7016125ad699bd7 100644 (file)
--- a/source/dofs/dof_tools.cc
+++ b/source/dofs/dof_tools.cc
@@ -687,6 +687,19 @@ namespace DoFTools
               {
                 const FiniteElement<dim, spacedim> &fe = cell->get_fe();
 
+                const auto reference_cell_type = cell->reference_cell_type();
+
+                const auto &cell_rc =
+                  ReferenceCell::internal::Info::get_cell(reference_cell_type);
+                const auto &face_rc =
+                  ReferenceCell::internal::Info::get_face(reference_cell_type,
+                                                          face);
+
+                const unsigned int n_vertices_per_cell = cell_rc.n_vertices();
+                const unsigned int n_lines_per_cell    = cell_rc.n_lines();
+                const unsigned int n_vertices_per_face = face_rc.n_vertices();
+                const unsigned int n_lines_per_face    = face_rc.n_lines();
+
                 const unsigned int dofs_per_face = fe.n_dofs_per_face(face);
                 face_dof_indices.resize(dofs_per_face);
                 cell->face(face)->get_dof_indices(face_dof_indices,
@@ -709,13 +722,26 @@ namespace DoFTools
                               (i < 2 * fe.n_dofs_per_vertex() ?
                                  i :
                                  i + 2 * fe.n_dofs_per_vertex()) :
-                              (dim == 3 ? (i < 4 * fe.n_dofs_per_vertex() ?
+                              (dim == 3 ? (i < n_vertices_per_face *
+                                                 fe.n_dofs_per_vertex() ?
                                              i :
-                                             (i < 4 * fe.n_dofs_per_vertex() +
-                                                    4 * fe.n_dofs_per_line() ?
-                                                i + 4 * fe.n_dofs_per_vertex() :
-                                                i + 4 * fe.n_dofs_per_vertex() +
-                                                  8 * fe.n_dofs_per_line())) :
+                                             (i < n_vertices_per_face *
+                                                      fe.n_dofs_per_vertex() +
+                                                    n_lines_per_face *
+                                                      fe.n_dofs_per_line() ?
+                                                (i - n_vertices_per_face *
+                                                       fe.n_dofs_per_vertex()) +
+                                                  n_vertices_per_cell *
+                                                    fe.n_dofs_per_vertex() :
+                                                (i -
+                                                 n_vertices_per_face *
+                                                   fe.n_dofs_per_vertex() -
+                                                 n_lines_per_face *
+                                                   fe.n_dofs_per_line()) +
+                                                  n_vertices_per_cell *
+                                                    fe.n_dofs_per_vertex() +
+                                                  n_lines_per_cell *
+                                                    fe.n_dofs_per_line())) :
                                           numbers::invalid_unsigned_int)));
                       if (fe.is_primitive(cell_index))
                         {

diff --git a/tests/simplex/extract_boundary_dofs.cc b/tests/simplex/extract_boundary_dofs.cc
new file mode 100644 (file)
index 0000000..711f944
--- /dev/null
+++ b/tests/simplex/extract_boundary_dofs.cc
@@ -0,0 +1,95 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2009 - 2018 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE.md at
+// the top level directory of deal.II.
+//
+// ---------------------------------------------------------------------
+
+
+
+// Like the test mpi/extract_boundary_dofs but for simplex meshes.
+
+#include <deal.II/distributed/shared_tria.h>
+
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/dofs/dof_tools.h>
+
+#include <deal.II/simplex/fe_lib.h>
+#include <deal.II/simplex/grid_generator.h>
+
+#include "../tests.h"
+
+
+
+template <int dim>
+void
+test()
+{
+  parallel::shared::Triangulation<dim> tr(
+    MPI_COMM_WORLD,
+    Triangulation<dim>::MeshSmoothing::none,
+    true,
+    parallel::shared::Triangulation<dim>::Settings::partition_zorder);
+
+  GridGenerator::subdivided_hyper_cube_with_simplices(tr, 4);
+
+  const Simplex::FE_P<dim> fe(2);
+  DoFHandler<dim>          dofh(tr);
+  dofh.distribute_dofs(fe);
+
+  IndexSet relevant_set, boundary_dofs;
+  DoFTools::extract_boundary_dofs(dofh,
+                                  std::vector<bool>(1, true),
+                                  boundary_dofs);
+  if (Utilities::MPI::this_mpi_process(MPI_COMM_WORLD) == 0)
+    boundary_dofs.write(deallog.get_file_stream());
+
+  // the result of extract_boundary_dofs is supposed to be a subset of the
+  // locally relevant dofs, so test this
+  DoFTools::extract_locally_relevant_dofs(dofh, relevant_set);
+  boundary_dofs.subtract_set(relevant_set);
+  AssertThrow(boundary_dofs.n_elements() == 0, ExcInternalError());
+}
+
+
+int
+main(int argc, char *argv[])
+{
+  Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1);
+
+  unsigned int myid = Utilities::MPI::this_mpi_process(MPI_COMM_WORLD);
+
+
+  deallog.push(Utilities::int_to_string(myid));
+
+  if (myid == 0)
+    {
+      initlog();
+
+      deallog.push("2d");
+      test<2>();
+      deallog.pop();
+
+      deallog.push("3d");
+      test<3>();
+      deallog.pop();
+    }
+  else
+    {
+      deallog.push("2d");
+      test<2>();
+      deallog.pop();
+
+      deallog.push("3d");
+      test<3>();
+      deallog.pop();
+    }
+}

diff --git a/tests/simplex/extract_boundary_dofs.mpirun=1.with_simplex_support=on.output b/tests/simplex/extract_boundary_dofs.mpirun=1.with_simplex_support=on.output
new file mode 100644 (file)
index 0000000..4bec103
--- /dev/null
+++ b/tests/simplex/extract_boundary_dofs.mpirun=1.with_simplex_support=on.output
@@ -0,0 +1,138 @@
+
+81 17
+0 4
+5 6
+9 11
+15 17
+21 23
+24 25
+26 28
+29 30
+42 43
+44 46
+47 48
+60 61
+62 64
+65 68
+70 72
+74 76
+78 81
+665 118
+0 11
+12 14
+17 20
+22 25
+26 34
+35 36
+37 38
+43 48
+49 51
+54 57
+60 70
+71 72
+73 76
+77 80
+81 86
+87 89
+94 96
+97 99
+105 108
+109 110
+116 121
+122 124
+125 126
+127 132
+133 135
+140 143
+144 147
+148 149
+155 157
+158 160
+166 171
+173 176
+177 180
+181 189
+190 191
+192 200
+201 202
+204 208
+209 211
+212 213
+214 224
+225 230
+231 233
+234 239
+244 245
+249 252
+255 258
+260 261
+266 268
+269 274
+275 276
+277 278
+284 287
+302 303
+305 308
+310 312
+313 315
+334 335
+336 338
+339 340
+341 344
+345 346
+348 352
+353 354
+355 356
+357 361
+362 363
+365 368
+369 377
+381 384
+386 389
+390 393
+395 396
+401 402
+406 409
+412 416
+417 418
+419 422
+424 426
+427 429
+448 449
+450 452
+453 454
+455 456
+462 465
+480 481
+483 486
+488 493
+494 495
+496 500
+501 502
+504 505
+506 507
+508 509
+510 516
+517 522
+523 525
+526 534
+535 538
+541 550
+552 556
+557 560
+561 571
+574 580
+582 585
+586 588
+591 595
+596 601
+602 604
+605 610
+611 613
+616 619
+621 624
+625 627
+628 636
+637 644
+645 653
+654 665

diff --git a/tests/simplex/extract_boundary_dofs.mpirun=5.with_simplex_support=on.output b/tests/simplex/extract_boundary_dofs.mpirun=5.with_simplex_support=on.output
new file mode 100644 (file)
index 0000000..3cd55a9
--- /dev/null
+++ b/tests/simplex/extract_boundary_dofs.mpirun=5.with_simplex_support=on.output
@@ -0,0 +1,102 @@
+
+81 8
+0 4
+5 6
+9 11
+15 17
+18 19
+20 21
+32 33
+34 37
+665 91
+0 11
+12 14
+17 20
+22 25
+26 34
+35 36
+37 38
+43 48
+49 51
+54 57
+60 70
+71 72
+73 76
+77 82
+83 85
+90 92
+93 95
+101 104
+105 106
+108 113
+114 116
+119 122
+123 124
+130 133
+134 141
+142 143
+144 152
+153 154
+156 160
+161 163
+164 169
+170 171
+175 178
+180 183
+185 186
+190 191
+192 193
+194 195
+196 197
+198 202
+203 204
+206 209
+210 211
+216 217
+223 226
+239 241
+242 244
+252 256
+260 263
+265 268
+269 272
+274 275
+277 279
+280 282
+304 306
+308 309
+310 313
+314 315
+372 373
+376 378
+380 384
+386 389
+414 417
+418 419
+421 422
+423 424
+425 431
+433 434
+473 476
+477 478
+486 487
+491 494
+497 501
+502 503
+579 582
+587 592
+593 595
+596 597
+598 600
+601 603
+609 614
+620 623
+624 626
+627 628
+629 631
+633 634
+641 642
+643 644
+645 648
+650 651
+660 664