]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Fix one testcase, see https://code.google.com/p/dealii/issues/detail?id=134 . Since...
authorWolfgang Bangerth <bangerth@math.tamu.edu>
Sun, 17 Nov 2013 19:35:33 +0000 (19:35 +0000)
committerWolfgang Bangerth <bangerth@math.tamu.edu>
Sun, 17 Nov 2013 19:35:33 +0000 (19:35 +0000)
git-svn-id: https://svn.dealii.org/trunk@31689 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/doc/news/changes.h
deal.II/include/deal.II/numerics/vector_tools.h
deal.II/include/deal.II/numerics/vector_tools.templates.h

index 38906d72b7c25f573990ac532d7439b07879fd43..04671e86902c9ca96f0e34d168c4ff9ab5a3b329 100644 (file)
@@ -94,6 +94,15 @@ inconvenience this causes.
 
 
 <ol>
+  <li>
+  Fixed: VectorTools::compute_no_normal_flux_constraints had a bug that
+  only appeared in rare cases at vertices of the domain if one adjacent
+  cell had two boundary indicators selected for no normal flux and another
+  had only one. This is now fixed.
+  <br>
+  (Wolfgang Bangerth, 2013/11/17)
+  </li>
+
   <li> New: introduced "make test" that runs a minimal set of tests. We
   encourage every user to run this, especially if they run in to problems.
   The tests are automatically picked depending on the configuration and
@@ -359,7 +368,7 @@ inconvenience this causes.
   New: All vector classes now have functions <code>extract_subvector_to()</code>
   that allow extracting not just a single value but a whole set.
   <br>
-  (Fahad Alrasched, 2013/09/02)
+  (Fahad Alrashed, 2013/09/02)
   </li>
 
   <li>
index d8c6be07ffc7e76ac02e5ec143636d3abb0aaf56..47e9823dac04128b4b5957321cfc8e042ce6ceec 100644 (file)
@@ -1404,7 +1404,7 @@ namespace VectorTools
 
 
   /**
-   * Compute the constraints that correspond to boundary conditions of the
+   * This function computes the constraints that correspond to boundary conditions of the
    * form $\vec n \cdot \vec u=0$, i.e. no normal flux if $\vec u$ is a
    * vector-valued quantity. These conditions have exactly the form handled by
    * the ConstraintMatrix class, so instead of creating a map between boundary
@@ -1563,8 +1563,19 @@ namespace VectorTools
    * There are cases where one cell contributes two tangential directions and
    * another one only one; for example, this would happen if both top and
    * front faces of the left cell belong to the boundary selected whereas only
-   * the top face of the right cell belongs to it. This case is not currently
-   * implemented.
+   * the top face of the right cell belongs to it, maybe indicating the the entire
+   * front part of the domain is a smooth manifold whereas the top really forms
+   * two separate manifolds that meet in a ridge, and that no-flux boundary
+   * conditions are only desired on the fron manifold and the right one on top.
+   * In cases like these, it's difficult to define what should happen. The
+   * current implementation simply ignores the one contribution from the
+   * cell that only contributes one normal vector. In the example shown, this
+   * is acceptable because the normal vector for the front face of the left
+   * cell is the same as the normal vector provided by the front face of
+   * the right cell (the surface is planar) but it would be a problem if the
+   * front manifold would be curved. Regardless, it is unclear how one would
+   * proceed in this case and ignoring the single cell is likely the best
+   * one can do.
    *
    *
    * <h4>Results</h4>
@@ -2032,7 +2043,7 @@ namespace VectorTools
    * point_difference() function.
    *
    * @note If the cell in which the point is found
-   * is not locally owned, an exception of type 
+   * is not locally owned, an exception of type
    * VectorTools::ExcPointNotAvailableHere
    * is thrown.
    */
@@ -2059,9 +2070,9 @@ namespace VectorTools
    * this function uses an
    * arbitrary mapping to evaluate
    * the difference.
-   * 
+   *
    * @note If the cell in which the point is found
-   * is not locally owned, an exception of type 
+   * is not locally owned, an exception of type
    * VectorTools::ExcPointNotAvailableHere
    * is thrown.
    */
@@ -2089,7 +2100,7 @@ namespace VectorTools
    * point_difference() function.
    *
    * @note If the cell in which the point is found
-   * is not locally owned, an exception of type 
+   * is not locally owned, an exception of type
    * VectorTools::ExcPointNotAvailableHere
    * is thrown.
    */
@@ -2104,7 +2115,7 @@ namespace VectorTools
   * Same as above for hp.
   *
   * @note If the cell in which the point is found
-  * is not locally owned, an exception of type 
+  * is not locally owned, an exception of type
   * VectorTools::ExcPointNotAvailableHere
   * is thrown.
   */
@@ -2134,7 +2145,7 @@ namespace VectorTools
    * "step-3".
    *
    * @note If the cell in which the point is found
-   * is not locally owned, an exception of type 
+   * is not locally owned, an exception of type
    * VectorTools::ExcPointNotAvailableHere
    * is thrown.
    */
@@ -2146,9 +2157,9 @@ namespace VectorTools
 
   /**
   * Same as above for hp.
-  * 
+  *
   * @note If the cell in which the point is found
-  * is not locally owned, an exception of type 
+  * is not locally owned, an exception of type
   * VectorTools::ExcPointNotAvailableHere
   * is thrown.
   */
@@ -2174,7 +2185,7 @@ namespace VectorTools
    * mapping to evaluate the difference.
    *
    * @note If the cell in which the point is found
-   * is not locally owned, an exception of type 
+   * is not locally owned, an exception of type
    * VectorTools::ExcPointNotAvailableHere
    * is thrown.
    */
@@ -2190,7 +2201,7 @@ namespace VectorTools
   * Same as above for hp.
   *
   * @note If the cell in which the point is found
-  * is not locally owned, an exception of type 
+  * is not locally owned, an exception of type
   * VectorTools::ExcPointNotAvailableHere
   * is thrown.
   */
@@ -2216,7 +2227,7 @@ namespace VectorTools
    * mapping to evaluate the difference.
    *
    * @note If the cell in which the point is found
-   * is not locally owned, an exception of type 
+   * is not locally owned, an exception of type
    * VectorTools::ExcPointNotAvailableHere
    * is thrown.
    */
@@ -2231,7 +2242,7 @@ namespace VectorTools
   * Same as above for hp.
   *
   * @note If the cell in which the point is found
-  * is not locally owned, an exception of type 
+  * is not locally owned, an exception of type
   * VectorTools::ExcPointNotAvailableHere
   * is thrown.
   */
index 7758845376fbcc9adad68c535296d100248f56f5..93e6847b74e0940da95d3c85f76a481706eb3ef8 100644 (file)
@@ -2404,6 +2404,16 @@ namespace VectorTools
     };
 
 
+    template <int dim>
+    std::ostream & operator << (std::ostream &out,
+                               const VectorDoFTuple<dim> &vdt)
+    {
+      for (unsigned int d=0; d<dim; ++d)
+       out << vdt.dof_indices[d] << (d < dim-1 ? " " : "");
+      return out;
+    }
+
+
 
     /**
      * Add the constraint
@@ -4382,6 +4392,12 @@ namespace VectorTools
                     .insert (std::make_pair (vector_dofs,
                                              std::make_pair (normal_vector,
                                                              cell)));
+#ifdef DEBUG_NO_NORMAL_FLUX
+                   std::cout << "Adding normal vector:" << std::endl
+                             << "   dofs=" << vector_dofs << std::endl
+                             << "   cell=" << cell << " at " << cell->center() << std::endl
+                             << "   normal=" << normal_vector << std::endl;
+#endif
                   }
             }
 
@@ -4409,6 +4425,18 @@ namespace VectorTools
         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;
+#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
@@ -4442,9 +4470,18 @@ namespace VectorTools
                 = 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;
+#endif
+
         // count the maximum number of contributions from each cell
         unsigned int max_n_contributions_per_cell = 1;
         for (typename CellToNormalsMap::const_iterator
@@ -4467,6 +4504,11 @@ namespace VectorTools
             // 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:
           {
 
@@ -4576,7 +4618,22 @@ namespace VectorTools
             // 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)
+            // 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
@@ -4586,7 +4643,22 @@ namespace VectorTools
                  contribution != cell_contributions.end();
                  ++contribution)
               {
-                Assert (contribution->second.size() == dim-1, ExcNotImplemented());
+#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
+
+               // as mentioned above, simply ignore cells that only
+               // contribute once
+                if (contribution->second.size() < dim-1)
+                 continue;
 
                 Tensor<1,dim> normals[dim-1];
                 {

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.