]> https://gitweb.dealii.org/ - dealii.git/commitdiff
simplex: face_interpolation for FE_Poly 11590/head
authorMarc Fehling <mafehling.git@gmail.com>
Fri, 15 Jan 2021 16:49:37 +0000 (09:49 -0700)
committerMarc Fehling <marc.fehling@gmx.net>
Thu, 28 Jan 2021 22:28:05 +0000 (15:28 -0700)
include/deal.II/simplex/fe_lib.h
source/base/qprojector.cc
source/fe/fe_q_base.cc
source/simplex/fe_lib.cc
tests/simplex/hanging_nodes_01.cc
tests/simplex/hanging_nodes_01.with_simplex_support=on.output
tests/simplex/hanging_nodes_02.cc
tests/simplex/hanging_nodes_02.with_simplex_support=on.output

index f8b8688a89dac7dd12d47f442b0b4d3fa38c4982..be7ce1cbb59b35c0ecb83e67643495af4d03503d 100644 (file)
@@ -62,6 +62,30 @@ namespace Simplex
       const RefinementCase<dim> &refinement_case =
         RefinementCase<dim>::isotropic_refinement) const override;
 
+    /**
+     * @copydoc dealii::FiniteElement::get_face_interpolation_matrix()
+     */
+    void
+    get_face_interpolation_matrix(const FiniteElement<dim, spacedim> &source_fe,
+                                  FullMatrix<double> &interpolation_matrix,
+                                  const unsigned int  face_no) const override;
+
+    /**
+     * @copydoc dealii::FiniteElement::get_subface_interpolation_matrix()
+     */
+    void
+    get_subface_interpolation_matrix(
+      const FiniteElement<dim, spacedim> &x_source_fe,
+      const unsigned int                  subface,
+      FullMatrix<double> &                interpolation_matrix,
+      const unsigned int                  face_no) const override;
+
+    /**
+     * @copydoc dealii::FiniteElement::hp_constraints_are_implemented()
+     */
+    bool
+    hp_constraints_are_implemented() const override;
+
     /**
      * @copydoc dealii::FiniteElement::convert_generalized_support_point_values_to_dof_values()
      */
index 6d05bf719b43ab4006d30820385bfcf616dfeef5..148223d578f83f6018fa2c96ac53e21d31fd7d00 100644 (file)
@@ -128,32 +128,57 @@ QProjector<2>::project_to_face(const ReferenceCell::Type reference_cell_type,
                                const unsigned int        face_no,
                                std::vector<Point<2>> &   q_points)
 {
-  Assert(reference_cell_type == ReferenceCell::Type::Quad, ExcNotImplemented());
-  (void)reference_cell_type;
-
   const unsigned int dim = 2;
   AssertIndexRange(face_no, GeometryInfo<dim>::faces_per_cell);
   Assert(q_points.size() == quadrature.size(),
          ExcDimensionMismatch(q_points.size(), quadrature.size()));
 
-  for (unsigned int p = 0; p < quadrature.size(); ++p)
-    switch (face_no)
-      {
-        case 0:
-          q_points[p] = Point<dim>(0, quadrature.point(p)(0));
-          break;
-        case 1:
-          q_points[p] = Point<dim>(1, quadrature.point(p)(0));
-          break;
-        case 2:
-          q_points[p] = Point<dim>(quadrature.point(p)(0), 0);
-          break;
-        case 3:
-          q_points[p] = Point<dim>(quadrature.point(p)(0), 1);
-          break;
-        default:
-          Assert(false, ExcInternalError());
-      }
+  if (reference_cell_type == ReferenceCell::Type::Tri)
+    {
+      // use linear polynomial to map the reference quadrature points correctly
+      // on faces, i.e., Simplex::ScalarPolynomial<1>(1)
+      for (unsigned int p = 0; p < quadrature.size(); ++p)
+        switch (face_no)
+          {
+            case 0:
+              q_points[p] = Point<dim>(quadrature.point(p)(0), 0);
+              break;
+            case 1:
+              q_points[p] =
+                Point<dim>(1 - quadrature.point(p)(0), quadrature.point(p)(0));
+              break;
+            case 2:
+              q_points[p] = Point<dim>(0, 1 - quadrature.point(p)(0));
+              break;
+            default:
+              Assert(false, ExcInternalError());
+          }
+    }
+  else if (reference_cell_type == ReferenceCell::Type::Quad)
+    {
+      for (unsigned int p = 0; p < quadrature.size(); ++p)
+        switch (face_no)
+          {
+            case 0:
+              q_points[p] = Point<dim>(0, quadrature.point(p)(0));
+              break;
+            case 1:
+              q_points[p] = Point<dim>(1, quadrature.point(p)(0));
+              break;
+            case 2:
+              q_points[p] = Point<dim>(quadrature.point(p)(0), 0);
+              break;
+            case 3:
+              q_points[p] = Point<dim>(quadrature.point(p)(0), 1);
+              break;
+            default:
+              Assert(false, ExcInternalError());
+          }
+    }
+  else
+    {
+      Assert(false, ExcInternalError());
+    }
 }
 
 
@@ -285,9 +310,6 @@ QProjector<2>::project_to_subface(const ReferenceCell::Type reference_cell_type,
                                   std::vector<Point<2>> &   q_points,
                                   const RefinementCase<1> &)
 {
-  Assert(reference_cell_type == ReferenceCell::Type::Quad, ExcNotImplemented());
-  (void)reference_cell_type;
-
   const unsigned int dim = 2;
   AssertIndexRange(face_no, GeometryInfo<dim>::faces_per_cell);
   AssertIndexRange(subface_no, GeometryInfo<dim>::max_children_per_face);
@@ -295,65 +317,130 @@ QProjector<2>::project_to_subface(const ReferenceCell::Type reference_cell_type,
   Assert(q_points.size() == quadrature.size(),
          ExcDimensionMismatch(q_points.size(), quadrature.size()));
 
-  for (unsigned int p = 0; p < quadrature.size(); ++p)
-    switch (face_no)
-      {
-        case 0:
-          switch (subface_no)
-            {
-              case 0:
-                q_points[p] = Point<dim>(0, quadrature.point(p)(0) / 2);
-                break;
-              case 1:
-                q_points[p] = Point<dim>(0, quadrature.point(p)(0) / 2 + 0.5);
-                break;
-              default:
-                Assert(false, ExcInternalError());
-            }
-          break;
-        case 1:
-          switch (subface_no)
-            {
-              case 0:
-                q_points[p] = Point<dim>(1, quadrature.point(p)(0) / 2);
-                break;
-              case 1:
-                q_points[p] = Point<dim>(1, quadrature.point(p)(0) / 2 + 0.5);
-                break;
-              default:
-                Assert(false, ExcInternalError());
-            }
-          break;
-        case 2:
-          switch (subface_no)
-            {
-              case 0:
-                q_points[p] = Point<dim>(quadrature.point(p)(0) / 2, 0);
-                break;
-              case 1:
-                q_points[p] = Point<dim>(quadrature.point(p)(0) / 2 + 0.5, 0);
-                break;
-              default:
-                Assert(false, ExcInternalError());
-            }
-          break;
-        case 3:
-          switch (subface_no)
-            {
-              case 0:
-                q_points[p] = Point<dim>(quadrature.point(p)(0) / 2, 1);
-                break;
-              case 1:
-                q_points[p] = Point<dim>(quadrature.point(p)(0) / 2 + 0.5, 1);
-                break;
-              default:
-                Assert(false, ExcInternalError());
-            }
-          break;
-
-        default:
-          Assert(false, ExcInternalError());
-      }
+  if (reference_cell_type == ReferenceCell::Type::Tri)
+    {
+      // use linear polynomial to map the reference quadrature points correctly
+      // on faces, i.e., Simplex::ScalarPolynomial<1>(1)
+      for (unsigned int p = 0; p < quadrature.size(); ++p)
+        switch (face_no)
+          {
+            case 0:
+              switch (subface_no)
+                {
+                  case 0:
+                    q_points[p] = Point<dim>(quadrature.point(p)(0) / 2, 0);
+                    break;
+                  case 1:
+                    q_points[p] =
+                      Point<dim>(0.5 + quadrature.point(p)(0) / 2, 0);
+                    break;
+                  default:
+                    Assert(false, ExcInternalError());
+                }
+              break;
+            case 1:
+              switch (subface_no)
+                {
+                  case 0:
+                    q_points[p] = Point<dim>(1 - quadrature.point(p)(0) / 2,
+                                             quadrature.point(p)(0) / 2);
+                    break;
+                  case 1:
+                    q_points[p] = Point<dim>(0.5 - quadrature.point(p)(0) / 2,
+                                             0.5 + quadrature.point(p)(0) / 2);
+                    break;
+                  default:
+                    Assert(false, ExcInternalError());
+                }
+              break;
+            case 2:
+              switch (subface_no)
+                {
+                  case 0:
+                    q_points[p] = Point<dim>(0, 1 - quadrature.point(p)(0) / 2);
+                    break;
+                  case 1:
+                    q_points[p] =
+                      Point<dim>(0, 0.5 - quadrature.point(p)(0) / 2);
+                    break;
+                  default:
+                    Assert(false, ExcInternalError());
+                }
+              break;
+            default:
+              Assert(false, ExcInternalError());
+          }
+    }
+  else if (reference_cell_type == ReferenceCell::Type::Quad)
+    {
+      for (unsigned int p = 0; p < quadrature.size(); ++p)
+        switch (face_no)
+          {
+            case 0:
+              switch (subface_no)
+                {
+                  case 0:
+                    q_points[p] = Point<dim>(0, quadrature.point(p)(0) / 2);
+                    break;
+                  case 1:
+                    q_points[p] =
+                      Point<dim>(0, quadrature.point(p)(0) / 2 + 0.5);
+                    break;
+                  default:
+                    Assert(false, ExcInternalError());
+                }
+              break;
+            case 1:
+              switch (subface_no)
+                {
+                  case 0:
+                    q_points[p] = Point<dim>(1, quadrature.point(p)(0) / 2);
+                    break;
+                  case 1:
+                    q_points[p] =
+                      Point<dim>(1, quadrature.point(p)(0) / 2 + 0.5);
+                    break;
+                  default:
+                    Assert(false, ExcInternalError());
+                }
+              break;
+            case 2:
+              switch (subface_no)
+                {
+                  case 0:
+                    q_points[p] = Point<dim>(quadrature.point(p)(0) / 2, 0);
+                    break;
+                  case 1:
+                    q_points[p] =
+                      Point<dim>(quadrature.point(p)(0) / 2 + 0.5, 0);
+                    break;
+                  default:
+                    Assert(false, ExcInternalError());
+                }
+              break;
+            case 3:
+              switch (subface_no)
+                {
+                  case 0:
+                    q_points[p] = Point<dim>(quadrature.point(p)(0) / 2, 1);
+                    break;
+                  case 1:
+                    q_points[p] =
+                      Point<dim>(quadrature.point(p)(0) / 2 + 0.5, 1);
+                    break;
+                  default:
+                    Assert(false, ExcInternalError());
+                }
+              break;
+
+            default:
+              Assert(false, ExcInternalError());
+          }
+    }
+  else
+    {
+      Assert(false, ExcInternalError());
+    }
 }
 
 
@@ -521,6 +608,17 @@ QProjector<2>::project_to_all_faces(
 {
   if (reference_cell_type == ReferenceCell::Type::Tri)
     {
+      const auto support_points_line =
+        [](const auto &face, const auto &orientation) -> std::vector<Point<2>> {
+        std::array<Point<2>, 2> vertices;
+        std::copy_n(face.first.begin(), face.first.size(), vertices.begin());
+        const auto temp =
+          ReferenceCell::Type(ReferenceCell::Type::Line)
+            .permute_according_orientation(vertices, orientation);
+        return std::vector<Point<2>>(temp.begin(),
+                                     temp.begin() + face.first.size());
+      };
+
       // reference faces (defined by its support points and arc length)
       const std::array<std::pair<std::array<Point<2>, 2>, double>, 3> faces = {
         {{{{Point<2>(0.0, 0.0), Point<2>(1.0, 0.0)}}, 1.0},
@@ -542,21 +640,10 @@ QProjector<2>::project_to_all_faces(
           {
             const auto &face = faces[face_no];
 
-            std::array<Point<2>, 2> support_points;
-
             // determine support point of the current line with the correct
             // orientation
-            switch (orientation)
-              {
-                case 0:
-                  support_points = {{face.first[1], face.first[0]}};
-                  break;
-                case 1:
-                  support_points = {{face.first[0], face.first[1]}};
-                  break;
-                default:
-                  Assert(false, ExcNotImplemented());
-              }
+            std::vector<Point<2>> support_points =
+              support_points_line(face, orientation);
 
             // the quadrature rule to be projected ...
             const auto &sub_quadrature_points =
index 77a219bca46b33cdf37ed38d54d874f00c5cd21f..ddbaed6fef8b40e9307a617f7df9470555767bd9 100644 (file)
@@ -697,6 +697,7 @@ FE_Q_Base<PolynomialType, dim, spacedim>::get_subface_interpolation_matrix(
             }
         }
 
+#ifdef DEBUG
       // make sure that the row sum of each of the matrices is 1 at this
       // point. this must be so since the shape functions sum up to 1
       for (unsigned int j = 0; j < source_fe->n_dofs_per_face(face_no); ++j)
@@ -708,6 +709,7 @@ FE_Q_Base<PolynomialType, dim, spacedim>::get_subface_interpolation_matrix(
 
           Assert(std::fabs(sum - 1) < eps, ExcInternalError());
         }
+#endif
     }
   else if (dynamic_cast<const FE_Nothing<dim> *>(&x_source_fe) != nullptr)
     {
index 61c14b5a43a880e0b5b0bd14d4874e2279d86806..2c31e505a81bec08e0a14a0ef47930ae85c381ee 100644 (file)
@@ -15,6 +15,8 @@
 
 #include <deal.II/base/config.h>
 
+#include <deal.II/base/qprojector.h>
+
 #include <deal.II/fe/fe_dgq.h>
 #include <deal.II/fe/fe_nothing.h>
 #include <deal.II/fe/fe_q.h>
@@ -370,6 +372,155 @@ namespace Simplex
 
 
 
+  template <int dim, int spacedim>
+  void
+  FE_Poly<dim, spacedim>::get_face_interpolation_matrix(
+    const FiniteElement<dim, spacedim> &x_source_fe,
+    FullMatrix<double> &                interpolation_matrix,
+    const unsigned int                  face_no) const
+  {
+    Assert(interpolation_matrix.m() == x_source_fe.n_dofs_per_face(face_no),
+           ExcDimensionMismatch(interpolation_matrix.m(),
+                                x_source_fe.n_dofs_per_face(face_no)));
+
+    if (const FE_Poly<dim, spacedim> *source_fe =
+          dynamic_cast<const FE_Poly<dim, spacedim> *>(&x_source_fe))
+      {
+        const Quadrature<dim - 1> quad_face_support(
+          source_fe->get_unit_face_support_points(face_no));
+
+        const double eps = 2e-13 * this->degree * (dim - 1);
+
+        std::vector<Point<dim>> face_quadrature_points(
+          quad_face_support.size());
+        QProjector<dim>::project_to_face(this->reference_cell_type(),
+                                         quad_face_support,
+                                         face_no,
+                                         face_quadrature_points);
+
+        for (unsigned int i = 0; i < source_fe->n_dofs_per_face(face_no); ++i)
+          for (unsigned int j = 0; j < this->n_dofs_per_face(face_no); ++j)
+            {
+              double matrix_entry =
+                this->shape_value(this->face_to_cell_index(j, 0),
+                                  face_quadrature_points[i]);
+
+              // Correct the interpolated value. I.e. if it is close to 1 or
+              // 0, make it exactly 1 or 0. Unfortunately, this is required to
+              // avoid problems with higher order elements.
+              if (std::fabs(matrix_entry - 1.0) < eps)
+                matrix_entry = 1.0;
+              if (std::fabs(matrix_entry) < eps)
+                matrix_entry = 0.0;
+
+              interpolation_matrix(i, j) = matrix_entry;
+            }
+
+#ifdef DEBUG
+        for (unsigned int j = 0; j < source_fe->n_dofs_per_face(face_no); ++j)
+          {
+            double sum = 0.;
+
+            for (unsigned int i = 0; i < this->n_dofs_per_face(face_no); ++i)
+              sum += interpolation_matrix(j, i);
+
+            Assert(std::fabs(sum - 1) < eps, ExcInternalError());
+          }
+#endif
+      }
+    else if (dynamic_cast<const FE_Nothing<dim> *>(&x_source_fe) != nullptr)
+      {
+        // nothing to do here, the FE_Nothing has no degrees of freedom anyway
+      }
+    else
+      AssertThrow(
+        false,
+        (typename FiniteElement<dim,
+                                spacedim>::ExcInterpolationNotImplemented()));
+  }
+
+
+
+  template <int dim, int spacedim>
+  void
+  FE_Poly<dim, spacedim>::get_subface_interpolation_matrix(
+    const FiniteElement<dim, spacedim> &x_source_fe,
+    const unsigned int                  subface,
+    FullMatrix<double> &                interpolation_matrix,
+    const unsigned int                  face_no) const
+  {
+    Assert(interpolation_matrix.m() == x_source_fe.n_dofs_per_face(face_no),
+           ExcDimensionMismatch(interpolation_matrix.m(),
+                                x_source_fe.n_dofs_per_face(face_no)));
+
+    if (const FE_Poly<dim, spacedim> *source_fe =
+          dynamic_cast<const FE_Poly<dim, spacedim> *>(&x_source_fe))
+      {
+        const Quadrature<dim - 1> quad_face_support(
+          source_fe->get_unit_face_support_points(face_no));
+
+        const double eps = 2e-13 * this->degree * (dim - 1);
+
+        std::vector<Point<dim>> subface_quadrature_points(
+          quad_face_support.size());
+        QProjector<dim>::project_to_subface(this->reference_cell_type(),
+                                            quad_face_support,
+                                            face_no,
+                                            subface,
+                                            subface_quadrature_points);
+
+        for (unsigned int i = 0; i < source_fe->n_dofs_per_face(face_no); ++i)
+          for (unsigned int j = 0; j < this->n_dofs_per_face(face_no); ++j)
+            {
+              double matrix_entry =
+                this->shape_value(this->face_to_cell_index(j, 0),
+                                  subface_quadrature_points[i]);
+
+              // Correct the interpolated value. I.e. if it is close to 1 or
+              // 0, make it exactly 1 or 0. Unfortunately, this is required to
+              // avoid problems with higher order elements.
+              if (std::fabs(matrix_entry - 1.0) < eps)
+                matrix_entry = 1.0;
+              if (std::fabs(matrix_entry) < eps)
+                matrix_entry = 0.0;
+
+              interpolation_matrix(i, j) = matrix_entry;
+            }
+
+#ifdef DEBUG
+        for (unsigned int j = 0; j < source_fe->n_dofs_per_face(face_no); ++j)
+          {
+            double sum = 0.;
+
+            for (unsigned int i = 0; i < this->n_dofs_per_face(face_no); ++i)
+              sum += interpolation_matrix(j, i);
+
+            Assert(std::fabs(sum - 1) < eps, ExcInternalError());
+          }
+#endif
+      }
+    else if (dynamic_cast<const FE_Nothing<dim> *>(&x_source_fe) != nullptr)
+      {
+        // nothing to do here, the FE_Nothing has no degrees of freedom anyway
+      }
+    else
+      AssertThrow(
+        false,
+        (typename FiniteElement<dim,
+                                spacedim>::ExcInterpolationNotImplemented()));
+  }
+
+
+
+  template <int dim, int spacedim>
+  bool
+  FE_Poly<dim, spacedim>::hp_constraints_are_implemented() const
+  {
+    return true;
+  }
+
+
+
   template <int dim, int spacedim>
   void
   FE_Poly<dim, spacedim>::
index 753ecdf45e1d21a94996c29e1726c92f1640b670..20762147c64157cdd1cdf387292e82108315d643 100644 (file)
 
 
 // verify hanging node constraints on locally h-refined simplex mesh
+//
+// dofs will be enumerated as follows
+//  scenario 1:
+//   1-------0
+//   |\      |
+//   |  \    |
+//   5---6   |
+//   |\  |\  |
+//   |  \|  \|
+//   3---4---2
+
 
 #include <deal.II/dofs/dof_handler.h>
 #include <deal.II/dofs/dof_tools.h>
 
+#include <deal.II/fe/mapping_fe.h>
+
 #include <deal.II/grid/grid_out.h>
 #include <deal.II/grid/tria.h>
 
 #include "../tests.h"
 
 
+// ----- diagnostics -----
+
+template <int dim>
+void
+print_dof_indices_on_faces(const DoFHandler<dim> &dofh)
+{
+  std::vector<types::global_dof_index> dof_indices;
+
+  for (const auto &cell : dofh.active_cell_iterators())
+    for (unsigned int f = 0; f < cell->n_faces(); ++f)
+      {
+        const auto &face = cell->face(f);
+
+        if (face->has_children())
+          {
+            for (unsigned int sf = 0; sf < face->n_children(); ++sf)
+              {
+                const auto &subface = face->child(sf);
+                Assert(subface->n_active_fe_indices() == 1, ExcInternalError());
+                const unsigned int subface_fe_index =
+                  subface->nth_active_fe_index(0);
+                const auto &subface_fe = dofh.get_fe(subface_fe_index);
+
+                dof_indices.resize(subface_fe.n_dofs_per_face(f));
+                subface->get_dof_indices(dof_indices, subface_fe_index);
+
+                deallog << "cell:" << cell->active_cell_index() << " face:" << f
+                        << " subface:" << sf << " dofs:";
+                for (const auto &i : dof_indices)
+                  deallog << i << " ";
+                deallog << std::endl;
+              }
+          }
+        else
+          {
+            Assert(face->n_active_fe_indices() == 1, ExcInternalError());
+            const unsigned int face_fe_index = face->nth_active_fe_index(0);
+            const auto &       face_fe       = dofh.get_fe(face_fe_index);
+
+            dof_indices.resize(face_fe.n_dofs_per_face(f));
+            face->get_dof_indices(dof_indices, face_fe_index);
+
+            deallog << "cell:" << cell->active_cell_index() << " face:" << f
+                    << " dofs:";
+            for (const auto &i : dof_indices)
+              deallog << i << " ";
+            deallog << std::endl;
+          }
+      }
+}
+
+
+template <int dim>
+void
+print_dof_points(const DoFHandler<dim> &dofh)
+{
+  std::vector<Point<dim>> points(dofh.n_dofs());
+  DoFTools::map_dofs_to_support_points(MappingFE<dim>(dofh.get_fe()),
+                                       dofh,
+                                       points);
+
+  for (unsigned int i = 0; i < dofh.n_dofs(); ++i)
+    deallog << "dof:" << i << " point:" << points[i] << std::endl;
+}
+
+
+
+// ----- test -----
+
 template <int dim>
 void
 test()
@@ -51,9 +133,14 @@ test()
   dofh.distribute_dofs(Simplex::FE_P<dim>(1));
   deallog << "ndofs: " << dofh.n_dofs() << std::endl;
 
+#if false
+  print_dof_points(dofh);
+  print_dof_indices_on_faces(dofh);
+#endif
+
   // hanging node constraints
   AffineConstraints<double> constraints;
-  // DoFTools::make_hanging_node_constraints(dofh, constraints);
+  DoFTools::make_hanging_node_constraints(dofh, constraints);
   constraints.print(deallog.get_file_stream());
 
   deallog << "OK" << std::endl;
index 1f032563725e69010bb3091f60d3d7556719cec1..8034a9e25676a3dfd446579df6af0d1bfc979325 100644 (file)
@@ -1,3 +1,5 @@
 
 DEAL:2d::ndofs: 7
+    6 2:  0.500000
+    6 1:  0.500000
 DEAL:2d::OK
index d3816830b14108dd194ab86a215fbfe378df70dc..024848424ed16f1f971be7270f9ee82ac8a80d58 100644 (file)
 
 
 // verify hanging node constraints on locally p-refined simplex mesh
+//
+// dofs will be enumerated as follows
+//  scenario 1:    scenario 2:
+//   6-------4      2---4---3
+//   |\      |      |\      |
+//   |  \    |      |  \    |
+//   3   2   |      |   5   6
+//   |    \  |      |    \  |
+//   |      \|      |      \|
+//   0---1---5      0-------1
+
 
 #include <deal.II/dofs/dof_handler.h>
 #include <deal.II/dofs/dof_tools.h>
 
+#include <deal.II/fe/mapping_fe.h>
+
 #include <deal.II/grid/tria.h>
 
 #include <deal.II/hp/fe_collection.h>
+#include <deal.II/hp/mapping_collection.h>
 
 #include <deal.II/lac/affine_constraints.h>
 
 #include "../tests.h"
 
 
+// ----- diagnostics -----
+
+template <int dim>
+void
+print_dof_indices_on_faces(const DoFHandler<dim> &dofh)
+{
+  std::vector<types::global_dof_index> dof_indices;
+
+  for (const auto &cell : dofh.active_cell_iterators())
+    for (unsigned int f = 0; f < cell->n_faces(); ++f)
+      {
+        const auto &face = cell->face(f);
+
+        Assert(!face->has_children(), ExcInternalError());
+
+        const unsigned int fe_index = cell->active_fe_index();
+        const auto &       fe       = cell->get_fe();
+
+        dof_indices.resize(fe.n_dofs_per_face(f));
+        face->get_dof_indices(dof_indices, fe_index);
+
+        deallog << "cell:" << cell->active_cell_index() << " face:" << f
+                << " dofs:";
+        for (const auto &i : dof_indices)
+          deallog << i << " ";
+        deallog << std::endl;
+      }
+}
+
+
+template <int dim>
+void
+print_dof_points(const DoFHandler<dim> &dofh)
+{
+  hp::MappingCollection<dim> mapping;
+  for (unsigned int i = 0; i < dofh.get_fe_collection().size(); ++i)
+    mapping.push_back(MappingFE<dim>(dofh.get_fe(i)));
+
+  std::vector<Point<dim>> points(dofh.n_dofs());
+  DoFTools::map_dofs_to_support_points(mapping, dofh, points);
+
+  for (unsigned int i = 0; i < dofh.n_dofs(); ++i)
+    deallog << "dof:" << i << " point:" << points[i] << std::endl;
+}
+
+
+// ----- test -----
+
 template <int dim>
 void
 test(const hp::FECollection<dim> &fes)
@@ -40,15 +102,25 @@ test(const hp::FECollection<dim> &fes)
   Triangulation<dim> tria;
   GridGenerator::subdivided_hyper_cube_with_simplices(tria, 1);
 
+#if false
+  GridOut grid_out;
+  grid_out.write_vtk(tria, deallog.get_file_stream());
+#endif
+
   DoFHandler<dim> dofh(tria);
   dofh.begin_active()->set_active_fe_index(1);
 
   dofh.distribute_dofs(fes);
   deallog << "ndofs: " << dofh.n_dofs() << std::endl;
 
+#if false
+  print_dof_points(dofh);
+  print_dof_indices_on_faces(dofh);
+#endif
+
   // hanging node constraints
   AffineConstraints<double> constraints;
-  // DoFTools::make_hanging_node_constraints(dofh, constraints);
+  DoFTools::make_hanging_node_constraints(dofh, constraints);
   constraints.print(deallog.get_file_stream());
 
   deallog << "OK" << std::endl;
index b905ef0c4e842b1573f2e4b696314d77f98a861c..0bc8c9bb95bfacfadc9f932e99c53528ef8346d7 100644 (file)
@@ -1,5 +1,9 @@
 
 DEAL:2d::ndofs: 7
+    2 6:  0.500000
+    2 5:  0.500000
 DEAL:2d::OK
 DEAL:2d::ndofs: 7
+    5 1:  0.500000
+    5 2:  0.500000
 DEAL:2d::OK

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.