]> https://gitweb.dealii.org/ - dealii.git/commitdiff
address comments 14167/head
authorJiaqi Zhang <jiaqi2@clemson.edu>
Wed, 24 Aug 2022 03:27:51 +0000 (23:27 -0400)
committerJiaqi Zhang <jiaqi2@clemson.edu>
Wed, 24 Aug 2022 03:27:51 +0000 (23:27 -0400)
include/deal.II/numerics/vector_tools_constraints.h
include/deal.II/numerics/vector_tools_constraints.templates.h
source/numerics/vector_tools_constraints.inst.in
tests/numerics/no_flux_17.cc
tests/numerics/no_flux_18.cc

index 9b2380d688e37347d9914d78448d408a4b43e0fb..8150a2b05d4c8b7e402b05d5fe62b886cbcde74a 100644 (file)
@@ -290,6 +290,33 @@ namespace VectorTools
          .template get_default_linear_mapping<dim, spacedim>()
 #else
          .ReferenceCell::get_default_linear_mapping<dim, spacedim>()
+#endif
+         ));
+
+  /**
+   * This function does the same as
+   * compute_nonzero_normal_flux_constraints(), but for the case of multigrid
+   * levels.
+   * @ingroup constraints
+   *
+   * @see
+   * @ref GlossBoundaryIndicator "Glossary entry on boundary indicators"
+   */
+  template <int dim, int spacedim>
+  void
+  compute_nonzero_normal_flux_constraints_on_level(
+    const DoFHandler<dim, spacedim> &   dof_handler,
+    const unsigned int                  first_vector_component,
+    const std::set<types::boundary_id> &boundary_ids,
+    const std::map<types::boundary_id, const Function<spacedim, double> *>
+      &                           function_map,
+    AffineConstraints<double> &   constraints,
+    const Mapping<dim, spacedim> &mapping =
+      (ReferenceCells::get_hypercube<dim>()
+#ifndef _MSC_VER
+         .template get_default_linear_mapping<dim, spacedim>()
+#else
+         .ReferenceCell::get_default_linear_mapping<dim, spacedim>()
 #endif
          ),
     const IndexSet &   refinement_edge_indices = IndexSet(),
@@ -323,6 +350,31 @@ namespace VectorTools
          .template get_default_linear_mapping<dim, spacedim>()
 #else
          .ReferenceCell::get_default_linear_mapping<dim, spacedim>()
+#endif
+         ));
+
+  /**
+   * This function does the same as
+   * compute_no_normal_flux_constraints(), but for the case of multigrid
+   * levels.
+   * @ingroup constraints
+   *
+   * @see
+   * @ref GlossBoundaryIndicator "Glossary entry on boundary indicators"
+   */
+  template <int dim, int spacedim>
+  void
+  compute_no_normal_flux_constraints_on_level(
+    const DoFHandler<dim, spacedim> &   dof_handler,
+    const unsigned int                  first_vector_component,
+    const std::set<types::boundary_id> &boundary_ids,
+    AffineConstraints<double> &         constraints,
+    const Mapping<dim, spacedim> &      mapping =
+      (ReferenceCells::get_hypercube<dim>()
+#ifndef _MSC_VER
+         .template get_default_linear_mapping<dim, spacedim>()
+#else
+         .ReferenceCell::get_default_linear_mapping<dim, spacedim>()
 #endif
          ),
     const IndexSet &   refinement_edge_indices = IndexSet(),
index e063d4ae628a20d928fd86dd7ebeaebe71073fe5..d15bd418176d502c71f34bfdd02c605ac4840104 100644 (file)
@@ -463,10 +463,13 @@ namespace VectorTools
     }
 
 
-
+    /**
+     * Compute the mappings from vector degrees of freedom to normal vectors @p dof_to_normals_map
+     * and vector degrees of freedom to prescribed normal fluxes @p dof_vector_to_b_values.
+     */
     template <int dim, int spacedim>
     void
-    get_dof_pairs(
+    map_dofs_to_normal_vectors_and_normal_fluxes(
       const typename DoFHandler<dim, spacedim>::cell_iterator &cell,
       const unsigned int                  first_vector_component,
       const std::set<types::boundary_id> &boundary_ids,
@@ -519,8 +522,8 @@ namespace VectorTools
                 // Refinement edge indices are going to be constrained to 0
                 // during a multigrid cycle and do not need no-normal-flux
                 // constraints, so skip them:
-                if (!refinement_edge_indices.is_element(face_dofs[i]) ||
-                    level == numbers::invalid_unsigned_int)
+                if (level == numbers::invalid_unsigned_int ||
+                    !refinement_edge_indices.is_element(face_dofs[i]))
                   {
                     // find corresponding other components of vector
                     internal::VectorDoFTuple<dim> vector_dofs;
@@ -634,485 +637,545 @@ namespace VectorTools
                   }
           }
     }
-  } // namespace internal
-
-
-
-  template <int dim, int spacedim>
-  void
-  compute_nonzero_normal_flux_constraints(
-    const DoFHandler<dim, spacedim> &   dof_handler,
-    const unsigned int                  first_vector_component,
-    const std::set<types::boundary_id> &boundary_ids,
-    const std::map<types::boundary_id, const Function<spacedim> *>
-      &                           function_map,
-    AffineConstraints<double> &   constraints,
-    const Mapping<dim, spacedim> &mapping,
-    const IndexSet &              refinement_edge_indices,
-    const unsigned int            level)
-  {
-    Assert(dim > 1,
-           ExcMessage("This function is not useful in 1d because it amounts "
-                      "to imposing Dirichlet values on the vector-valued "
-                      "quantity."));
-
-    // create FE and mapping collections for all elements in use by this
-    // DoFHandler
-    const hp::FECollection<dim, spacedim> &fe_collection =
-      dof_handler.get_fe_collection();
-    hp::MappingCollection<dim, spacedim> mapping_collection;
-    for (unsigned int i = 0; i < fe_collection.size(); ++i)
-      mapping_collection.push_back(mapping);
 
-    // TODO: the implementation makes the assumption that all faces have the
-    // same number of dofs
-    AssertDimension(dof_handler.get_fe().n_unique_faces(), 1);
-    const unsigned int face_no = 0;
 
-    // now also create a quadrature collection for the faces of a cell. fill
-    // it with a quadrature formula with the support points on faces for each
-    // FE
-    hp::QCollection<dim - 1> face_quadrature_collection;
-    for (unsigned int i = 0; i < fe_collection.size(); ++i)
-      {
-        const std::vector<Point<dim - 1>> &unit_support_points =
-          fe_collection[i].get_unit_face_support_points(face_no);
 
-        Assert(unit_support_points.size() ==
-                 fe_collection[i].n_dofs_per_face(face_no),
-               ExcInternalError());
+    /**
+     * This is the internal function that computes the nonzero normal
+     * flux constraints on active cells
+     * if @p level is an invalid unsigned integer or level cells if the cell level is provided.
+     * It's called by compute_nonzero_normal_flux_constraints() and
+     * compute_nonzero_normal_flux_constraints_on_level() so as to have
+     * separate interfaces for the active and level cells.
+     */
+    template <int dim, int spacedim>
+    void
+    compute_nonzero_normal_flux_constraints_active_or_level(
+      const DoFHandler<dim, spacedim> &   dof_handler,
+      const unsigned int                  first_vector_component,
+      const std::set<types::boundary_id> &boundary_ids,
+      const std::map<types::boundary_id, const Function<spacedim> *>
+        &                           function_map,
+      AffineConstraints<double> &   constraints,
+      const Mapping<dim, spacedim> &mapping,
+      const IndexSet &              refinement_edge_indices = IndexSet(),
+      const unsigned int            level = numbers::invalid_unsigned_int)
+    {
+      Assert(dim > 1,
+             ExcMessage("This function is not useful in 1d because it amounts "
+                        "to imposing Dirichlet values on the vector-valued "
+                        "quantity."));
+
+      // create FE and mapping collections for all elements in use by this
+      // DoFHandler
+      const hp::FECollection<dim, spacedim> &fe_collection =
+        dof_handler.get_fe_collection();
+      hp::MappingCollection<dim, spacedim> mapping_collection;
+      for (unsigned int i = 0; i < fe_collection.size(); ++i)
+        mapping_collection.push_back(mapping);
+
+      // TODO: the implementation makes the assumption that all faces have the
+      // same number of dofs
+      AssertDimension(dof_handler.get_fe().n_unique_faces(), 1);
+      const unsigned int face_no = 0;
+
+      // now also create a quadrature collection for the faces of a cell. fill
+      // it with a quadrature formula with the support points on faces for each
+      // FE
+      hp::QCollection<dim - 1> face_quadrature_collection;
+      for (unsigned int i = 0; i < fe_collection.size(); ++i)
+        {
+          const std::vector<Point<dim - 1>> &unit_support_points =
+            fe_collection[i].get_unit_face_support_points(face_no);
 
-        face_quadrature_collection.push_back(
-          Quadrature<dim - 1>(unit_support_points));
-      }
+          Assert(unit_support_points.size() ==
+                   fe_collection[i].n_dofs_per_face(face_no),
+                 ExcInternalError());
 
-    // now create the object with which we will generate the normal vectors
-    hp::FEFaceValues<dim, spacedim> x_fe_face_values(mapping_collection,
-                                                     fe_collection,
-                                                     face_quadrature_collection,
-                                                     update_quadrature_points |
-                                                       update_normal_vectors);
+          face_quadrature_collection.push_back(
+            Quadrature<dim - 1>(unit_support_points));
+        }
 
-    // have a map that stores normal vectors for each vector-dof tuple we want
-    // to constrain. since we can get at the same vector dof tuple more than
-    // once (for example if it is located at a vertex that we visit from all
-    // adjacent cells), we will want to average later on the normal vectors
-    // computed on different cells as described in the documentation of this
-    // function. however, we can only average if the contributions came from
-    // different cells, whereas we want to constrain twice or more in case the
-    // contributions came from different faces of the same cell
-    // (i.e. constrain not just the *average normal direction* but *all normal
-    // directions* we find). consequently, we also have to store which cell a
-    // normal vector was computed on
-    using DoFToNormalsMap = std::multimap<
-      internal::VectorDoFTuple<dim>,
-      std::pair<Tensor<1, dim>,
-                typename DoFHandler<dim, spacedim>::cell_iterator>>;
-    std::map<internal::VectorDoFTuple<dim>, Vector<double>>
-      dof_vector_to_b_values;
+      // now create the object with which we will generate the normal vectors
+      hp::FEFaceValues<dim, spacedim> x_fe_face_values(
+        mapping_collection,
+        fe_collection,
+        face_quadrature_collection,
+        update_quadrature_points | update_normal_vectors);
+
+      // have a map that stores normal vectors for each vector-dof tuple we want
+      // to constrain. since we can get at the same vector dof tuple more than
+      // once (for example if it is located at a vertex that we visit from all
+      // adjacent cells), we will want to average later on the normal vectors
+      // computed on different cells as described in the documentation of this
+      // function. however, we can only average if the contributions came from
+      // different cells, whereas we want to constrain twice or more in case the
+      // contributions came from different faces of the same cell
+      // (i.e. constrain not just the *average normal direction* but *all normal
+      // directions* we find). consequently, we also have to store which cell a
+      // normal vector was computed on
+      using DoFToNormalsMap = std::multimap<
+        internal::VectorDoFTuple<dim>,
+        std::pair<Tensor<1, dim>,
+                  typename DoFHandler<dim, spacedim>::cell_iterator>>;
+      std::map<internal::VectorDoFTuple<dim>, Vector<double>>
+        dof_vector_to_b_values;
 
-    DoFToNormalsMap dof_to_normals_map;
+      DoFToNormalsMap dof_to_normals_map;
 
-    const unsigned int n_dof = dof_handler.n_dofs();
+      const unsigned int n_dof = dof_handler.n_dofs();
 
-    if (level == numbers::invalid_unsigned_int)
-      {
-        // active cells
-        for (const auto &cell : dof_handler.active_cell_iterators())
-          if (!cell->is_artificial())
-            {
-              internal::get_dof_pairs(cell,
-                                      first_vector_component,
-                                      boundary_ids,
-                                      function_map,
-                                      x_fe_face_values,
-                                      n_dof,
-                                      refinement_edge_indices,
-                                      level,
-                                      dof_to_normals_map,
-                                      dof_vector_to_b_values);
-            }
-      }
-    else
-      {
-        // level cells
-        for (const auto &cell : dof_handler.cell_iterators_on_level(level))
-          if (cell->level_subdomain_id() != numbers::artificial_subdomain_id &&
-              cell->level_subdomain_id() != numbers::invalid_subdomain_id)
-            {
-              internal::get_dof_pairs(cell,
-                                      first_vector_component,
-                                      boundary_ids,
-                                      function_map,
-                                      x_fe_face_values,
-                                      n_dof,
-                                      refinement_edge_indices,
-                                      level,
-                                      dof_to_normals_map,
-                                      dof_vector_to_b_values);
-            }
-      }
+      if (level == numbers::invalid_unsigned_int)
+        {
+          // active cells
+          for (const auto &cell : dof_handler.active_cell_iterators())
+            if (!cell->is_artificial())
+              {
+                internal::map_dofs_to_normal_vectors_and_normal_fluxes(
+                  cell,
+                  first_vector_component,
+                  boundary_ids,
+                  function_map,
+                  x_fe_face_values,
+                  n_dof,
+                  refinement_edge_indices,
+                  level,
+                  dof_to_normals_map,
+                  dof_vector_to_b_values);
+              }
+        }
+      else
+        {
+          // level cells
+          for (const auto &cell : dof_handler.cell_iterators_on_level(level))
+            if (cell->level_subdomain_id() !=
+                  numbers::artificial_subdomain_id &&
+                cell->level_subdomain_id() != numbers::invalid_subdomain_id)
+              {
+                internal::map_dofs_to_normal_vectors_and_normal_fluxes(
+                  cell,
+                  first_vector_component,
+                  boundary_ids,
+                  function_map,
+                  x_fe_face_values,
+                  n_dof,
+                  refinement_edge_indices,
+                  level,
+                  dof_to_normals_map,
+                  dof_vector_to_b_values);
+              }
+        }
 
-    // Now do something with the collected information. To this end, loop
-    // through all sets of pairs (dofs,normal_vector) and identify which
-    // entries belong to the same set of dofs and then do as described in the
-    // documentation, i.e. either average the normal vector or don't for this
-    // particular set of dofs
-    typename DoFToNormalsMap::const_iterator p = dof_to_normals_map.begin();
+      // Now do something with the collected information. To this end, loop
+      // through all sets of pairs (dofs,normal_vector) and identify which
+      // entries belong to the same set of dofs and then do as described in the
+      // documentation, i.e. either average the normal vector or don't for this
+      // particular set of dofs
+      typename DoFToNormalsMap::const_iterator p = dof_to_normals_map.begin();
 
-    while (p != dof_to_normals_map.end())
-      {
-        // first find the range of entries in the multimap that corresponds to
-        // the same vector-dof tuple. as usual, we define the range
-        // half-open. the first entry of course is 'p'
-        typename DoFToNormalsMap::const_iterator same_dof_range[2] = {p};
-        for (++p; p != dof_to_normals_map.end(); ++p)
-          if (p->first != same_dof_range[0]->first)
-            {
-              same_dof_range[1] = p;
-              break;
-            }
-        if (p == dof_to_normals_map.end())
-          same_dof_range[1] = dof_to_normals_map.end();
+      while (p != dof_to_normals_map.end())
+        {
+          // first find the range of entries in the multimap that corresponds to
+          // the same vector-dof tuple. as usual, we define the range
+          // half-open. the first entry of course is 'p'
+          typename DoFToNormalsMap::const_iterator same_dof_range[2] = {p};
+          for (++p; p != dof_to_normals_map.end(); ++p)
+            if (p->first != same_dof_range[0]->first)
+              {
+                same_dof_range[1] = p;
+                break;
+              }
+          if (p == dof_to_normals_map.end())
+            same_dof_range[1] = dof_to_normals_map.end();
 
 #ifdef DEBUG_NO_NORMAL_FLUX
-        std::cout << "For dof indices <" << p->first
-                  << ">, found the following normals" << std::endl;
-        for (typename DoFToNormalsMap::const_iterator q = same_dof_range[0];
-             q != same_dof_range[1];
-             ++q)
-          std::cout << "   " << q->second.first << " from cell "
-                    << q->second.second << std::endl;
+          std::cout << "For dof indices <" << p->first
+                    << ">, found the following normals" << std::endl;
+          for (typename DoFToNormalsMap::const_iterator q = same_dof_range[0];
+               q != same_dof_range[1];
+               ++q)
+            std::cout << "   " << q->second.first << " from cell "
+                      << q->second.second << std::endl;
 #endif
 
 
-        // now compute the reverse mapping: for each of the cells that
-        // contributed to the current set of vector dofs, add up the normal
-        // vectors. the values of the map are pairs of normal vectors and
-        // number of cells that have contributed
-        using CellToNormalsMap =
-          std::map<typename DoFHandler<dim, spacedim>::cell_iterator,
-                   std::pair<Tensor<1, dim>, unsigned int>>;
-
-        CellToNormalsMap cell_to_normals_map;
-        for (typename DoFToNormalsMap::const_iterator q = same_dof_range[0];
-             q != same_dof_range[1];
-             ++q)
-          if (cell_to_normals_map.find(q->second.second) ==
-              cell_to_normals_map.end())
-            cell_to_normals_map[q->second.second] =
-              std::make_pair(q->second.first, 1U);
-          else
-            {
-              const Tensor<1, dim> old_normal =
-                cell_to_normals_map[q->second.second].first;
-              const unsigned int old_count =
-                cell_to_normals_map[q->second.second].second;
-
-              Assert(old_count > 0, ExcInternalError());
-
-              // in the same entry, store again the now averaged normal vector
-              // and the new count
+          // now compute the reverse mapping: for each of the cells that
+          // contributed to the current set of vector dofs, add up the normal
+          // vectors. the values of the map are pairs of normal vectors and
+          // number of cells that have contributed
+          using CellToNormalsMap =
+            std::map<typename DoFHandler<dim, spacedim>::cell_iterator,
+                     std::pair<Tensor<1, dim>, unsigned int>>;
+
+          CellToNormalsMap cell_to_normals_map;
+          for (typename DoFToNormalsMap::const_iterator q = same_dof_range[0];
+               q != same_dof_range[1];
+               ++q)
+            if (cell_to_normals_map.find(q->second.second) ==
+                cell_to_normals_map.end())
               cell_to_normals_map[q->second.second] =
-                std::make_pair((old_normal * old_count + q->second.first) /
-                                 (old_count + 1),
-                               old_count + 1);
-            }
-        Assert(cell_to_normals_map.size() >= 1, ExcInternalError());
+                std::make_pair(q->second.first, 1U);
+            else
+              {
+                const Tensor<1, dim> old_normal =
+                  cell_to_normals_map[q->second.second].first;
+                const unsigned int old_count =
+                  cell_to_normals_map[q->second.second].second;
+
+                Assert(old_count > 0, ExcInternalError());
+
+                // in the same entry, store again the now averaged normal vector
+                // and the new count
+                cell_to_normals_map[q->second.second] =
+                  std::make_pair((old_normal * old_count + q->second.first) /
+                                   (old_count + 1),
+                                 old_count + 1);
+              }
+          Assert(cell_to_normals_map.size() >= 1, ExcInternalError());
 
 #ifdef DEBUG_NO_NORMAL_FLUX
-        std::cout << "   cell_to_normals_map:" << std::endl;
-        for (typename CellToNormalsMap::const_iterator x =
-               cell_to_normals_map.begin();
-             x != cell_to_normals_map.end();
-             ++x)
-          std::cout << "      " << x->first << " -> (" << x->second.first << ','
-                    << x->second.second << ')' << std::endl;
+          std::cout << "   cell_to_normals_map:" << std::endl;
+          for (typename CellToNormalsMap::const_iterator x =
+                 cell_to_normals_map.begin();
+               x != cell_to_normals_map.end();
+               ++x)
+            std::cout << "      " << x->first << " -> (" << x->second.first
+                      << ',' << x->second.second << ')' << std::endl;
 #endif
 
-        // count the maximum number of contributions from each cell
-        unsigned int max_n_contributions_per_cell = 1;
-        for (typename CellToNormalsMap::const_iterator x =
-               cell_to_normals_map.begin();
-             x != cell_to_normals_map.end();
-             ++x)
-          max_n_contributions_per_cell =
-            std::max(max_n_contributions_per_cell, x->second.second);
-
-        // verify that each cell can have only contributed at most dim times,
-        // since that is the maximum number of faces that come together at a
-        // single place
-        Assert(max_n_contributions_per_cell <= dim, ExcInternalError());
-
-        switch (max_n_contributions_per_cell)
-          {
-            // first deal with the case that a number of cells all have
-            // registered that they have a normal vector defined at the
-            // location of a given vector dof, and that each of them have
-            // encountered this vector dof exactly once while looping over all
-            // their faces. as stated in the documentation, this is the case
-            // where we want to simply average over all normal vectors
-            //
-            // the typical case is in 2d where multiple cells meet at one
-            // vertex sitting on the boundary. same in 3d for a vertex that
-            // is associated with only one of the boundary indicators passed
-            // to this function
-            case 1:
-              {
-                // compute the average normal vector from all the ones that
-                // have the same set of dofs. we could add them up and divide
-                // them by the number of additions, or simply normalize them
-                // right away since we want them to have unit length anyway
-                Tensor<1, dim> normal;
-                for (typename CellToNormalsMap::const_iterator x =
-                       cell_to_normals_map.begin();
-                     x != cell_to_normals_map.end();
-                     ++x)
-                  normal += x->second.first;
-                normal /= normal.norm();
-
-                // normalize again
-                for (unsigned int d = 0; d < dim; ++d)
-                  if (std::fabs(normal[d]) < 1e-13)
-                    normal[d] = 0;
-                normal /= normal.norm();
-
-                // then construct constraints from this:
-                const internal::VectorDoFTuple<dim> &dof_indices =
-                  same_dof_range[0]->first;
-                double               normal_value = 0.;
-                const Vector<double> b_values =
-                  dof_vector_to_b_values[dof_indices];
-                for (unsigned int i = 0; i < dim; ++i)
-                  normal_value += b_values[i] * normal[i];
-                internal::add_constraint(dof_indices,
-                                         normal,
-                                         constraints,
-                                         normal_value);
-
-                break;
-              }
-
-            // this is the slightly more complicated case that a single cell
-            // has contributed with exactly DIM normal vectors to the same set
-            // of vector dofs. this is what happens in a corner in 2d and 3d
-            // (but not on an edge in 3d, where we have only 2, i.e. <DIM,
-            // contributions. Here we do not want to average the normal
-            // vectors. Since we have DIM contributions, let's assume (and
-            // verify) that they are in fact all linearly independent; in that
-            // case, all vector components are constrained and we need to set
-            // all of them to the corresponding boundary values
-            case dim:
-              {
-                // assert that indeed only a single cell has contributed
-                Assert(cell_to_normals_map.size() == 1, ExcInternalError());
-
-                // check linear independence by computing the determinant of
-                // the matrix created from all the normal vectors. if they are
-                // linearly independent, then the determinant is nonzero. if
-                // they are orthogonal, then the matrix is in fact equal to 1
-                // (since they are all unit vectors); make sure the
-                // determinant is larger than 1e-3 to avoid cases where cells
-                // are degenerate
+          // count the maximum number of contributions from each cell
+          unsigned int max_n_contributions_per_cell = 1;
+          for (typename CellToNormalsMap::const_iterator x =
+                 cell_to_normals_map.begin();
+               x != cell_to_normals_map.end();
+               ++x)
+            max_n_contributions_per_cell =
+              std::max(max_n_contributions_per_cell, x->second.second);
+
+          // verify that each cell can have only contributed at most dim times,
+          // since that is the maximum number of faces that come together at a
+          // single place
+          Assert(max_n_contributions_per_cell <= dim, ExcInternalError());
+
+          switch (max_n_contributions_per_cell)
+            {
+              // first deal with the case that a number of cells all have
+              // registered that they have a normal vector defined at the
+              // location of a given vector dof, and that each of them have
+              // encountered this vector dof exactly once while looping over all
+              // their faces. as stated in the documentation, this is the case
+              // where we want to simply average over all normal vectors
+              //
+              // the typical case is in 2d where multiple cells meet at one
+              // vertex sitting on the boundary. same in 3d for a vertex that
+              // is associated with only one of the boundary indicators passed
+              // to this function
+              case 1:
                 {
-                  Tensor<2, dim> t;
-
-                  typename DoFToNormalsMap::const_iterator x =
-                    same_dof_range[0];
-                  for (unsigned int i = 0; i < dim; ++i, ++x)
-                    for (unsigned int j = 0; j < dim; ++j)
-                      t[i][j] = x->second.first[j];
-
-                  Assert(
-                    std::fabs(determinant(t)) > 1e-3,
-                    ExcMessage(
-                      "Found a set of normal vectors that are nearly collinear."));
+                  // compute the average normal vector from all the ones that
+                  // have the same set of dofs. we could add them up and divide
+                  // them by the number of additions, or simply normalize them
+                  // right away since we want them to have unit length anyway
+                  Tensor<1, dim> normal;
+                  for (typename CellToNormalsMap::const_iterator x =
+                         cell_to_normals_map.begin();
+                       x != cell_to_normals_map.end();
+                       ++x)
+                    normal += x->second.first;
+                  normal /= normal.norm();
+
+                  // normalize again
+                  for (unsigned int d = 0; d < dim; ++d)
+                    if (std::fabs(normal[d]) < 1e-13)
+                      normal[d] = 0;
+                  normal /= normal.norm();
+
+                  // then construct constraints from this:
+                  const internal::VectorDoFTuple<dim> &dof_indices =
+                    same_dof_range[0]->first;
+                  double               normal_value = 0.;
+                  const Vector<double> b_values =
+                    dof_vector_to_b_values[dof_indices];
+                  for (unsigned int i = 0; i < dim; ++i)
+                    normal_value += b_values[i] * normal[i];
+                  internal::add_constraint(dof_indices,
+                                           normal,
+                                           constraints,
+                                           normal_value);
+
+                  break;
                 }
 
-                // so all components of this vector dof are constrained. enter
-                // this into the AffineConstraints object
-                //
-                // ignore dofs already constrained
-                const internal::VectorDoFTuple<dim> &dof_indices =
-                  same_dof_range[0]->first;
-                const Vector<double> b_values =
-                  dof_vector_to_b_values[dof_indices];
-                for (unsigned int i = 0; i < dim; ++i)
-                  if (!constraints.is_constrained(
-                        same_dof_range[0]->first.dof_indices[i]) &&
-                      constraints.can_store_line(
-                        same_dof_range[0]->first.dof_indices[i]))
-                    {
-                      const types::global_dof_index line =
-                        dof_indices.dof_indices[i];
-                      constraints.add_line(line);
-                      if (std::fabs(b_values[i]) >
-                          std::numeric_limits<double>::epsilon())
-                        constraints.set_inhomogeneity(line, b_values[i]);
-                      // no add_entries here
-                    }
+              // this is the slightly more complicated case that a single cell
+              // has contributed with exactly DIM normal vectors to the same set
+              // of vector dofs. this is what happens in a corner in 2d and 3d
+              // (but not on an edge in 3d, where we have only 2, i.e. <DIM,
+              // contributions. Here we do not want to average the normal
+              // vectors. Since we have DIM contributions, let's assume (and
+              // verify) that they are in fact all linearly independent; in that
+              // case, all vector components are constrained and we need to set
+              // all of them to the corresponding boundary values
+              case dim:
+                {
+                  // assert that indeed only a single cell has contributed
+                  Assert(cell_to_normals_map.size() == 1, ExcInternalError());
+
+                  // check linear independence by computing the determinant of
+                  // the matrix created from all the normal vectors. if they are
+                  // linearly independent, then the determinant is nonzero. if
+                  // they are orthogonal, then the matrix is in fact equal to 1
+                  // (since they are all unit vectors); make sure the
+                  // determinant is larger than 1e-3 to avoid cases where cells
+                  // are degenerate
+                  {
+                    Tensor<2, dim> t;
 
-                break;
-              }
+                    typename DoFToNormalsMap::const_iterator x =
+                      same_dof_range[0];
+                    for (unsigned int i = 0; i < dim; ++i, ++x)
+                      for (unsigned int j = 0; j < dim; ++j)
+                        t[i][j] = x->second.first[j];
 
-            // this is the case of an edge contribution in 3d, i.e. the vector
-            // is constrained in two directions but not the third.
-            default:
-              {
-                Assert(dim >= 3, ExcNotImplemented());
-                Assert(max_n_contributions_per_cell == 2, ExcInternalError());
-
-                // as described in the documentation, let us first collect
-                // what each of the cells contributed at the current point. we
-                // use a std::list instead of a std::set (which would be more
-                // natural) because std::set requires that the stored elements
-                // are comparable with operator<
-                using CellContributions =
-                  std::map<typename DoFHandler<dim, spacedim>::cell_iterator,
-                           std::list<Tensor<1, dim>>>;
-                CellContributions cell_contributions;
-
-                for (typename DoFToNormalsMap::const_iterator q =
-                       same_dof_range[0];
-                     q != same_dof_range[1];
-                     ++q)
-                  cell_contributions[q->second.second].push_back(
-                    q->second.first);
-                Assert(cell_contributions.size() >= 1, ExcInternalError());
-
-                // now for each cell that has contributed determine the number
-                // of normal vectors it has contributed. we currently only
-                // implement if this is dim-1 for all cells (if a single cell
-                // has contributed dim, or if all adjacent cells have
-                // contributed 1 normal vector, this is already handled
-                // above).
-                //
-                // we only implement the case that all cells contribute
-                // dim-1 because we assume that we are following an edge
-                // of the domain (think: we are looking at a vertex
-                // located on one of the edges of a refined cube where the
-                // boundary indicators of the two adjacent faces of the
-                // cube are both listed in the set of boundary indicators
-                // passed to this function). in that case, all cells along
-                // that edge of the domain are assumed to have contributed
-                // dim-1 normal vectors. however, there are cases where
-                // this assumption is not justified (see the lengthy
-                // explanation in test no_flux_12.cc) and in those cases
-                // we simply ignore the cell that contributes only
-                // once. this is also discussed at length in the
-                // documentation of this function.
-                //
-                // for each contributing cell compute the tangential vector
-                // that remains unconstrained
-                std::list<Tensor<1, dim>> tangential_vectors;
-                for (typename CellContributions::const_iterator contribution =
-                       cell_contributions.begin();
-                     contribution != cell_contributions.end();
-                     ++contribution)
-                  {
-#ifdef DEBUG_NO_NORMAL_FLUX
-                    std::cout
-                      << "   Treating edge case with dim-1 contributions."
-                      << std::endl
-                      << "   Looking at cell " << contribution->first
-                      << " which has contributed these normal vectors:"
-                      << std::endl;
-                    for (typename std::list<Tensor<1, dim>>::const_iterator t =
-                           contribution->second.begin();
-                         t != contribution->second.end();
-                         ++t)
-                      std::cout << "      " << *t << std::endl;
-#endif
+                    Assert(
+                      std::fabs(determinant(t)) > 1e-3,
+                      ExcMessage(
+                        "Found a set of normal vectors that are nearly collinear."));
+                  }
 
-                    // as mentioned above, simply ignore cells that only
-                    // contribute once
-                    if (contribution->second.size() < dim - 1)
-                      continue;
+                  // so all components of this vector dof are constrained. enter
+                  // this into the AffineConstraints object
+                  //
+                  // ignore dofs already constrained
+                  const internal::VectorDoFTuple<dim> &dof_indices =
+                    same_dof_range[0]->first;
+                  const Vector<double> b_values =
+                    dof_vector_to_b_values[dof_indices];
+                  for (unsigned int i = 0; i < dim; ++i)
+                    if (!constraints.is_constrained(
+                          same_dof_range[0]->first.dof_indices[i]) &&
+                        constraints.can_store_line(
+                          same_dof_range[0]->first.dof_indices[i]))
+                      {
+                        const types::global_dof_index line =
+                          dof_indices.dof_indices[i];
+                        constraints.add_line(line);
+                        if (std::fabs(b_values[i]) >
+                            std::numeric_limits<double>::epsilon())
+                          constraints.set_inhomogeneity(line, b_values[i]);
+                        // no add_entries here
+                      }
+
+                  break;
+                }
 
-                    Tensor<1, dim> normals[dim - 1];
+              // this is the case of an edge contribution in 3d, i.e. the vector
+              // is constrained in two directions but not the third.
+              default:
+                {
+                  Assert(dim >= 3, ExcNotImplemented());
+                  Assert(max_n_contributions_per_cell == 2, ExcInternalError());
+
+                  // as described in the documentation, let us first collect
+                  // what each of the cells contributed at the current point. we
+                  // use a std::list instead of a std::set (which would be more
+                  // natural) because std::set requires that the stored elements
+                  // are comparable with operator<
+                  using CellContributions =
+                    std::map<typename DoFHandler<dim, spacedim>::cell_iterator,
+                             std::list<Tensor<1, dim>>>;
+                  CellContributions cell_contributions;
+
+                  for (typename DoFToNormalsMap::const_iterator q =
+                         same_dof_range[0];
+                       q != same_dof_range[1];
+                       ++q)
+                    cell_contributions[q->second.second].push_back(
+                      q->second.first);
+                  Assert(cell_contributions.size() >= 1, ExcInternalError());
+
+                  // now for each cell that has contributed determine the number
+                  // of normal vectors it has contributed. we currently only
+                  // implement if this is dim-1 for all cells (if a single cell
+                  // has contributed dim, or if all adjacent cells have
+                  // contributed 1 normal vector, this is already handled
+                  // above).
+                  //
+                  // we only implement the case that all cells contribute
+                  // dim-1 because we assume that we are following an edge
+                  // of the domain (think: we are looking at a vertex
+                  // located on one of the edges of a refined cube where the
+                  // boundary indicators of the two adjacent faces of the
+                  // cube are both listed in the set of boundary indicators
+                  // passed to this function). in that case, all cells along
+                  // that edge of the domain are assumed to have contributed
+                  // dim-1 normal vectors. however, there are cases where
+                  // this assumption is not justified (see the lengthy
+                  // explanation in test no_flux_12.cc) and in those cases
+                  // we simply ignore the cell that contributes only
+                  // once. this is also discussed at length in the
+                  // documentation of this function.
+                  //
+                  // for each contributing cell compute the tangential vector
+                  // that remains unconstrained
+                  std::list<Tensor<1, dim>> tangential_vectors;
+                  for (typename CellContributions::const_iterator contribution =
+                         cell_contributions.begin();
+                       contribution != cell_contributions.end();
+                       ++contribution)
                     {
-                      unsigned int index = 0;
+#ifdef DEBUG_NO_NORMAL_FLUX
+                      std::cout
+                        << "   Treating edge case with dim-1 contributions."
+                        << std::endl
+                        << "   Looking at cell " << contribution->first
+                        << " which has contributed these normal vectors:"
+                        << std::endl;
                       for (typename std::list<Tensor<1, dim>>::const_iterator
                              t = contribution->second.begin();
                            t != contribution->second.end();
-                           ++t, ++index)
-                        normals[index] = *t;
-                      Assert(index == dim - 1, ExcInternalError());
-                    }
+                           ++t)
+                        std::cout << "      " << *t << std::endl;
+#endif
+
+                      // as mentioned above, simply ignore cells that only
+                      // contribute once
+                      if (contribution->second.size() < dim - 1)
+                        continue;
 
-                    // calculate the tangent as the outer product of the
-                    // normal vectors. since these vectors do not need to be
-                    // orthogonal (think, for example, the case of the
-                    // deal.II/no_flux_07 test: a sheared cube in 3d, with Q2
-                    // elements, where we have constraints from the two normal
-                    // vectors of two faces of the sheared cube that are not
-                    // perpendicular to each other), we have to normalize the
-                    // outer product
-                    Tensor<1, dim> tangent;
-                    switch (dim)
+                      Tensor<1, dim> normals[dim - 1];
                       {
-                        case 3:
-                          // take cross product between normals[0] and
-                          // normals[1]. write it in the current form (with
-                          // [dim-2]) to make sure that compilers don't warn
-                          // about out-of-bounds accesses -- the warnings are
-                          // bogus since we get here only for dim==3, but at
-                          // least one isn't quite smart enough to notice this
-                          // and warns when compiling the function in 2d
-                          tangent =
-                            cross_product_3d(normals[0], normals[dim - 2]);
-                          break;
-                        default:
-                          Assert(false, ExcNotImplemented());
+                        unsigned int index = 0;
+                        for (typename std::list<Tensor<1, dim>>::const_iterator
+                               t = contribution->second.begin();
+                             t != contribution->second.end();
+                             ++t, ++index)
+                          normals[index] = *t;
+                        Assert(index == dim - 1, ExcInternalError());
                       }
 
-                    Assert(
-                      std::fabs(tangent.norm()) > 1e-12,
-                      ExcMessage(
-                        "Two normal vectors from adjacent faces are almost "
-                        "parallel."));
-                    tangent /= tangent.norm();
+                      // calculate the tangent as the outer product of the
+                      // normal vectors. since these vectors do not need to be
+                      // orthogonal (think, for example, the case of the
+                      // deal.II/no_flux_07 test: a sheared cube in 3d, with Q2
+                      // elements, where we have constraints from the two normal
+                      // vectors of two faces of the sheared cube that are not
+                      // perpendicular to each other), we have to normalize the
+                      // outer product
+                      Tensor<1, dim> tangent;
+                      switch (dim)
+                        {
+                          case 3:
+                            // take cross product between normals[0] and
+                            // normals[1]. write it in the current form (with
+                            // [dim-2]) to make sure that compilers don't warn
+                            // about out-of-bounds accesses -- the warnings are
+                            // bogus since we get here only for dim==3, but at
+                            // least one isn't quite smart enough to notice this
+                            // and warns when compiling the function in 2d
+                            tangent =
+                              cross_product_3d(normals[0], normals[dim - 2]);
+                            break;
+                          default:
+                            Assert(false, ExcNotImplemented());
+                        }
+
+                      Assert(
+                        std::fabs(tangent.norm()) > 1e-12,
+                        ExcMessage(
+                          "Two normal vectors from adjacent faces are almost "
+                          "parallel."));
+                      tangent /= tangent.norm();
+
+                      tangential_vectors.push_back(tangent);
+                    }
 
-                    tangential_vectors.push_back(tangent);
+                  // go through the list of tangents and make sure that they all
+                  // roughly point in the same direction as the first one (i.e.
+                  // have an angle less than 90 degrees); if they don't then
+                  // flip their sign
+                  {
+                    const Tensor<1, dim> first_tangent =
+                      tangential_vectors.front();
+                    typename std::list<Tensor<1, dim>>::iterator t =
+                      tangential_vectors.begin();
+                    ++t;
+                    for (; t != tangential_vectors.end(); ++t)
+                      if (*t * first_tangent < 0)
+                        *t *= -1;
                   }
 
-                // go through the list of tangents and make sure that they all
-                // roughly point in the same direction as the first one (i.e.
-                // have an angle less than 90 degrees); if they don't then
-                // flip their sign
-                {
-                  const Tensor<1, dim> first_tangent =
-                    tangential_vectors.front();
-                  typename std::list<Tensor<1, dim>>::iterator t =
-                    tangential_vectors.begin();
-                  ++t;
-                  for (; t != tangential_vectors.end(); ++t)
-                    if (*t * first_tangent < 0)
-                      *t *= -1;
+                  // now compute the average tangent and normalize it
+                  Tensor<1, dim> average_tangent;
+                  for (typename std::list<Tensor<1, dim>>::const_iterator t =
+                         tangential_vectors.begin();
+                       t != tangential_vectors.end();
+                       ++t)
+                    average_tangent += *t;
+                  average_tangent /= average_tangent.norm();
+
+                  // now all that is left is that we add the constraints that
+                  // the vector is parallel to the tangent
+                  const internal::VectorDoFTuple<dim> &dof_indices =
+                    same_dof_range[0]->first;
+                  const Vector<double> b_values =
+                    dof_vector_to_b_values[dof_indices];
+                  internal::add_tangentiality_constraints(dof_indices,
+                                                          average_tangent,
+                                                          constraints,
+                                                          b_values);
                 }
+            }
+        }
+    }
 
-                // now compute the average tangent and normalize it
-                Tensor<1, dim> average_tangent;
-                for (typename std::list<Tensor<1, dim>>::const_iterator t =
-                       tangential_vectors.begin();
-                     t != tangential_vectors.end();
-                     ++t)
-                  average_tangent += *t;
-                average_tangent /= average_tangent.norm();
-
-                // now all that is left is that we add the constraints that
-                // the vector is parallel to the tangent
-                const internal::VectorDoFTuple<dim> &dof_indices =
-                  same_dof_range[0]->first;
-                const Vector<double> b_values =
-                  dof_vector_to_b_values[dof_indices];
-                internal::add_tangentiality_constraints(dof_indices,
-                                                        average_tangent,
-                                                        constraints,
-                                                        b_values);
-              }
-          }
-      }
+  } // namespace internal
+
+
+
+  template <int dim, int spacedim>
+  void
+  compute_nonzero_normal_flux_constraints(
+    const DoFHandler<dim, spacedim> &   dof_handler,
+    const unsigned int                  first_vector_component,
+    const std::set<types::boundary_id> &boundary_ids,
+    const std::map<types::boundary_id, const Function<spacedim> *>
+      &                           function_map,
+    AffineConstraints<double> &   constraints,
+    const Mapping<dim, spacedim> &mapping)
+  {
+    internal::compute_nonzero_normal_flux_constraints_active_or_level(
+      dof_handler,
+      first_vector_component,
+      boundary_ids,
+      function_map,
+      constraints,
+      mapping);
+  }
+
+
+
+  template <int dim, int spacedim>
+  void
+  compute_nonzero_normal_flux_constraints_on_level(
+    const DoFHandler<dim, spacedim> &   dof_handler,
+    const unsigned int                  first_vector_component,
+    const std::set<types::boundary_id> &boundary_ids,
+    const std::map<types::boundary_id, const Function<spacedim> *>
+      &                           function_map,
+    AffineConstraints<double> &   constraints,
+    const Mapping<dim, spacedim> &mapping,
+    const IndexSet &              refinement_edge_indices,
+    const unsigned int            level)
+  {
+    internal::compute_nonzero_normal_flux_constraints_active_or_level(
+      dof_handler,
+      first_vector_component,
+      boundary_ids,
+      function_map,
+      constraints,
+      mapping,
+      refinement_edge_indices,
+      level);
   }
 
   namespace internal
@@ -1350,6 +1413,30 @@ namespace VectorTools
   template <int dim, int spacedim>
   void
   compute_no_normal_flux_constraints(
+    const DoFHandler<dim, spacedim> &   dof_handler,
+    const unsigned int                  first_vector_component,
+    const std::set<types::boundary_id> &boundary_ids,
+    AffineConstraints<double> &         constraints,
+    const Mapping<dim, spacedim> &      mapping)
+  {
+    Functions::ZeroFunction<dim>                             zero_function(dim);
+    std::map<types::boundary_id, const Function<spacedim> *> function_map;
+    for (const types::boundary_id boundary_id : boundary_ids)
+      function_map[boundary_id] = &zero_function;
+    internal::compute_nonzero_normal_flux_constraints_active_or_level(
+      dof_handler,
+      first_vector_component,
+      boundary_ids,
+      function_map,
+      constraints,
+      mapping);
+  }
+
+
+
+  template <int dim, int spacedim>
+  void
+  compute_no_normal_flux_constraints_on_level(
     const DoFHandler<dim, spacedim> &   dof_handler,
     const unsigned int                  first_vector_component,
     const std::set<types::boundary_id> &boundary_ids,
@@ -1362,14 +1449,15 @@ namespace VectorTools
     std::map<types::boundary_id, const Function<spacedim> *> function_map;
     for (const types::boundary_id boundary_id : boundary_ids)
       function_map[boundary_id] = &zero_function;
-    compute_nonzero_normal_flux_constraints(dof_handler,
-                                            first_vector_component,
-                                            boundary_ids,
-                                            function_map,
-                                            constraints,
-                                            mapping,
-                                            refinement_edge_indices,
-                                            level);
+    internal::compute_nonzero_normal_flux_constraints_active_or_level(
+      dof_handler,
+      first_vector_component,
+      boundary_ids,
+      function_map,
+      constraints,
+      mapping,
+      refinement_edge_indices,
+      level);
   }
 
 
index 9e47b54f84b94afa8ed9c0869cf16dd18076a061..2f659f92f17fd41e7a55f16bd0c58ebb9dab5fba 100644 (file)
@@ -26,6 +26,16 @@ for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS)
 #  if deal_II_dimension != 1
       template void
       compute_nonzero_normal_flux_constraints(
+        const DoFHandler<deal_II_dimension> &dof_handler,
+        const unsigned int                   first_vector_component,
+        const std::set<types::boundary_id> & boundary_ids,
+        const std::map<types::boundary_id, const Function<deal_II_dimension> *>
+          &                               function_map,
+        AffineConstraints<double> &       constraints,
+        const Mapping<deal_II_dimension> &mapping);
+
+      template void
+      compute_nonzero_normal_flux_constraints_on_level(
         const DoFHandler<deal_II_dimension> &dof_handler,
         const unsigned int                   first_vector_component,
         const std::set<types::boundary_id> & boundary_ids,
@@ -49,6 +59,14 @@ for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS)
 
       template void
       compute_no_normal_flux_constraints(
+        const DoFHandler<deal_II_dimension> &dof_handler,
+        const unsigned int                   first_vector_component,
+        const std::set<types::boundary_id> & boundary_ids,
+        AffineConstraints<double> &          constraints,
+        const Mapping<deal_II_dimension> &   mapping);
+
+      template void
+      compute_no_normal_flux_constraints_on_level(
         const DoFHandler<deal_II_dimension> &dof_handler,
         const unsigned int                   first_vector_component,
         const std::set<types::boundary_id> & boundary_ids,
index edd36d3a766183a801cd96a5219bb5e07980ca50..7536dfddc8761038879efe334939f65416e2cb97 100644 (file)
@@ -15,7 +15,7 @@
 
 
 
-// add tests for compute_no_normal_flux_constraints_on_level()
+// test compute_no_normal_flux_constraints_on_lefel().
 
 #include <deal.II/base/exceptions.h>
 #include <deal.II/base/function.h>
@@ -172,13 +172,14 @@ get_constraints_on_levels(
       const IndexSet &refinement_edge_indices =
         mg_constrained_dofs.get_refinement_edge_indices(level);
 
-      VectorTools::compute_no_normal_flux_constraints(dof_handler,
-                                                      0,
-                                                      no_flux_boundary,
-                                                      user_level_constraints,
-                                                      mapping,
-                                                      refinement_edge_indices,
-                                                      level);
+      VectorTools::compute_no_normal_flux_constraints_on_level(
+        dof_handler,
+        0,
+        no_flux_boundary,
+        user_level_constraints,
+        mapping,
+        refinement_edge_indices,
+        level);
       user_level_constraints.close();
       level_constraints[level].copy_from(user_level_constraints);
       // deallog<<" ***level cell constraints:
index 6813fdec13682cfe461bf60ea5a30a9c8c8d7d1c..5bd134fdcfbbd0da9498664d071d331fb59db337 100644 (file)
@@ -1133,13 +1133,14 @@ namespace StokesClass
         const auto &  refinement_edge_indices =
           mg_constrained_dofs.get_refinement_edge_indices(level);
 
-        VectorTools::compute_no_normal_flux_constraints(dof_handler_u,
-                                                        0,
-                                                        dirichlet_boundary,
-                                                        level_constraints,
-                                                        mapping,
-                                                        refinement_edge_indices,
-                                                        level);
+        VectorTools::compute_no_normal_flux_constraints_on_level(
+          dof_handler_u,
+          0,
+          dirichlet_boundary,
+          level_constraints,
+          mapping,
+          refinement_edge_indices,
+          level);
         level_constraints.close();
 
         mg_constrained_dofs.add_user_constraints(level, level_constraints);

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.