]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Implement multiway DoF identities on vertices/lines/quads.
authorWolfgang Bangerth <bangerth@colostate.edu>
Fri, 20 Aug 2021 19:18:22 +0000 (13:18 -0600)
committerWolfgang Bangerth <bangerth@colostate.edu>
Tue, 24 Aug 2021 19:25:21 +0000 (13:25 -0600)
include/deal.II/hp/fe_collection.h
source/hp/fe_collection.cc

index 57f73ebd184f9317ee8ccc076ec72ef2eda83d5e..a5d4ba63260fb83f9e1c65df41513d49d424ceb6 100644 (file)
@@ -300,6 +300,38 @@ namespace hp
     bool
     hp_constraints_are_implemented() const;
 
+    /**
+     * This function combines the functionality of the
+     * FiniteElement::hp_vertex_dof_identities() into multi-way comparisons.
+     * Given a set of elements (whose indices are provided as argument), this
+     * function determines identities between degrees of freedom of these
+     * elements at a vertex.
+     *
+     * The function returns a vector of such identities, where each element of
+     * the vector is a set of pairs `(fe_index,dof_index)` that identifies
+     * the `fe_index` (an element of the `fes` argument to this function) of
+     * an element and the `dof_index` indicates the how-manyth degree of freedom
+     * of that element on a vertex participates in this identity.
+     */
+    std::vector<std::set<std::pair<unsigned int, unsigned int>>>
+    hp_vertex_dof_identities(const std::set<unsigned int> &fes) const;
+
+    /**
+     * Same as hp_vertex_dof_indices(), except that the function treats degrees
+     * of freedom on lines.
+     */
+    std::vector<std::set<std::pair<unsigned int, unsigned int>>>
+    hp_line_dof_identities(const std::set<unsigned int> &fes) const;
+
+    /**
+     * Same as hp_vertex_dof_indices(), except that the function treats degrees
+     * of freedom on quads.
+     */
+    std::vector<std::set<std::pair<unsigned int, unsigned int>>>
+    hp_quad_dof_identities(const std::set<unsigned int> &fes,
+                           const unsigned int            face_no = 0) const;
+
+
     /**
      * Return the indices of finite elements in this FECollection that dominate
      * all elements associated with the provided set of indices @p fes.
index 80661cd0891da7147a3d9cd717fde06e11873e25..8b08beac8573f3807b2d742ba1ec9f64c574880f 100644 (file)
 
 #include <deal.II/hp/fe_collection.h>
 
+#include <set>
+
+
+
 DEAL_II_NAMESPACE_OPEN
 
 namespace hp
@@ -286,6 +290,206 @@ namespace hp
 
 
 
+  namespace
+  {
+    /**
+     * Implement the action of the hp_*_dof_identities() functions
+     * in a generic way.
+     */
+    std::vector<std::set<std::pair<unsigned int, unsigned int>>>
+    compute_hp_dof_identities(
+      const std::set<unsigned int> &fes,
+      const std::function<std::vector<std::pair<unsigned int, unsigned int>>(
+        const unsigned int,
+        const unsigned int)> &      query_identities)
+    {
+      // Consider all degrees of freedom of the identified elements (represented
+      // via (fe_index,dof_index) pairs) as the nodes in a graph. Then draw
+      // edges for all DoFs that are identified based on what the elements
+      // selected in the argument say. Let us first build this graph, where we
+      // only store the edges of the graph, and as a consequence ignore nodes
+      // (DoFs) that simply don't show up at all in any of the identities:
+      using Node  = std::pair<unsigned int, unsigned int>;
+      using Edge  = std::pair<Node, Node>;
+      using Graph = std::set<Edge>;
+
+      Graph identities_graph;
+      for (const unsigned int fe_index_1 : fes)
+        for (const unsigned int fe_index_2 : fes)
+          if (fe_index_1 != fe_index_2)
+            for (const auto &identity :
+                 query_identities(fe_index_1, fe_index_2))
+              identities_graph.emplace(Edge(Node(fe_index_1, identity.first),
+                                            Node(fe_index_2, identity.second)));
+
+#ifdef DEBUG
+      // Now verify that indeed the graph is symmetric: If one element
+      // declares that certain ones of its DoFs are to be unified with those
+      // of the other, then the other one should agree with this. As a
+      // consequence of this test succeeding, we know that the graph is actually
+      // undirected.
+      for (const auto &edge : identities_graph)
+        Assert(identities_graph.find({edge.second, edge.first}) !=
+                 identities_graph.end(),
+               ExcInternalError());
+#endif
+
+      // The next step is that we ought to verify that if there is an identity
+      // between (fe1,dof1) and (fe2,dof2), as well as with (fe2,dof2) and
+      // (fe3,dof3), then there should also be an identity between (fe1,dof1)
+      // and (fe3,dof3). The same logic should apply to chains of four
+      // identities.
+      //
+      // This means that the graph we have built above must be composed of a
+      // collection of complete sub-graphs (complete = each possible edge in the
+      // sub-graph exists) -- or, using a different term, that the graph
+      // consists of a number of "cliques". Each of these cliques is then one
+      // extended identity between two or more DoFs, and these are the ones that
+      // we will want to return.
+      //
+      // To ascertain that this is true, and to figure out what we want to
+      // return, we need to divide the graph into its sub-graphs and then ensure
+      // that each sub-graph is indeed a clique. This is slightly cumbersome,
+      // but can be done as follows:
+      // - pick one edge 'e' of G
+      // - add e=(n1,n2) to the sub-graph SG
+      // - set N={n1,n2}
+      // - loop over the remainder of the edges 'e' of the graph:
+      //   - if 'e' has one or both nodes in N:
+      //     - add 'e' to SG and
+      //     - add its two nodes to N (they may already be in there)
+      //     - remove 'e' from G
+      //
+      // In general, this may not find the entire sub-graph if the edges are
+      // stored in random order. For example, if the graph consisted of the
+      // following edges in this order:
+      //   (a,b)
+      //   (c,d)
+      //   (a,c)
+      //   (a,d)
+      //   (b,c)
+      //   (b,d)
+      // then the graph itself is a clique, but the algorithm outlined above
+      // would skip the edge (c,d) because neither of the nodes are already
+      // in the set N which at that point is still (a,b).
+      //
+      // But, we store the graph with a std::set, which stored edges in sorted
+      // order where the order is the lexicographic order of nodes. This ensures
+      // that we really capture all edges that correspond to a sub-graph (but
+      // we will assert this as well).
+      //
+      // (For programmatic reasons, we skip the removal of 'e' from G in a first
+      // run through because it modifies the graph and thus invalidates
+      // iterators. But because SG stores all of these edges, we can remove them
+      // all from G after collecting the edges in SG.)
+      std::vector<std::set<Node>> identities;
+      while (identities_graph.size() > 0)
+        {
+          Graph          sub_graph;       // SG
+          std::set<Node> sub_graph_nodes; // N
+
+          sub_graph.emplace(*identities_graph.begin());
+          sub_graph_nodes.emplace(identities_graph.begin()->first);
+          sub_graph_nodes.emplace(identities_graph.begin()->second);
+
+          for (const Edge &e : identities_graph)
+            if ((sub_graph_nodes.find(e.first) != sub_graph_nodes.end()) ||
+                (sub_graph_nodes.find(e.second) != sub_graph_nodes.end()))
+              {
+                sub_graph.insert(e);
+                sub_graph_nodes.insert(e.first);
+                sub_graph_nodes.insert(e.second);
+              }
+
+          // We have now obtained a sub-graph from the overall graph.
+          // Now remove it from the bigger graph
+          for (const Edge &e : sub_graph)
+            identities_graph.erase(e);
+
+#if DEBUG
+          // There are three checks we ought to perform:
+          // - That the sub-graph is undirected, i.e. that every edge appears
+          //   in both directions
+          for (const auto &edge : sub_graph)
+            Assert(sub_graph.find({edge.second, edge.first}) != sub_graph.end(),
+                   ExcInternalError());
+
+          // - None of the nodes in the sub-graph should have appeared in
+          //   any of the other sub-graphs. If they did, then we have a bug
+          //   in extracting sub-graphs. This is actually more easily checked
+          //   the other way around: none of the nodes of the sub-graph we
+          //   just extracted should be in any of the edges of the *remaining*
+          //   graph
+          for (const Node &n : sub_graph_nodes)
+            for (const Edge &e : identities_graph)
+              Assert((n != e.first) && (n != e.second), ExcInternalError());
+          // - Second, the sub-graph we just extracted needs to be complete,
+          //   i.e.,
+          //   be a "clique". We check this by counting how many edges it has.
+          //   for 'n' nodes in 'N', we need to have n*(n-1) edges (we store
+          //   both directed edges).
+          Assert(sub_graph.size() ==
+                   sub_graph_nodes.size() * (sub_graph_nodes.size() - 1),
+                 ExcInternalError());
+#endif
+
+          // At this point we're sure that we have extracted a complete
+          // sub-graph ("clique"). The DoFs involved are all identical then, and
+          // we will store this identity so we can return it later.
+          identities.emplace_back(std::move(sub_graph_nodes));
+        }
+
+      return identities;
+    }
+  } // namespace
+
+
+
+  template <int dim, int spacedim>
+  std::vector<std::set<std::pair<unsigned int, unsigned int>>>
+  FECollection<dim, spacedim>::hp_vertex_dof_identities(
+    const std::set<unsigned int> &fes) const
+  {
+    auto query_vertex_dof_identities = [this](const unsigned int fe_index_1,
+                                              const unsigned int fe_index_2) {
+      return (*this)[fe_index_1].hp_vertex_dof_identities((*this)[fe_index_2]);
+    };
+    return compute_hp_dof_identities(fes, query_vertex_dof_identities);
+  }
+
+
+
+  template <int dim, int spacedim>
+  std::vector<std::set<std::pair<unsigned int, unsigned int>>>
+  FECollection<dim, spacedim>::hp_line_dof_identities(
+    const std::set<unsigned int> &fes) const
+  {
+    auto query_line_dof_identities = [this](const unsigned int fe_index_1,
+                                            const unsigned int fe_index_2) {
+      return (*this)[fe_index_1].hp_line_dof_identities((*this)[fe_index_2]);
+    };
+    return compute_hp_dof_identities(fes, query_line_dof_identities);
+  }
+
+
+
+  template <int dim, int spacedim>
+  std::vector<std::set<std::pair<unsigned int, unsigned int>>>
+  FECollection<dim, spacedim>::hp_quad_dof_identities(
+    const std::set<unsigned int> &fes,
+    const unsigned int            face_no) const
+  {
+    auto query_quad_dof_identities = [this,
+                                      face_no](const unsigned int fe_index_1,
+                                               const unsigned int fe_index_2) {
+      return (*this)[fe_index_1].hp_quad_dof_identities((*this)[fe_index_2],
+                                                        face_no);
+    };
+    return compute_hp_dof_identities(fes, query_quad_dof_identities);
+  }
+
+
+
   template <int dim, int spacedim>
   void
   FECollection<dim, spacedim>::set_hierarchy(

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.