From 0277432238e8d7ea73e60d96b32579bb5eeccd03 Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Sun, 17 Nov 2013 19:35:33 +0000 Subject: [PATCH] Fix one testcase, see https://code.google.com/p/dealii/issues/detail?id=134 . Since we are having so much trouble with this function, I also left my debug code in there, guarded with \#ifdef DEBUG_NO_NORMAL_FLUX. git-svn-id: https://svn.dealii.org/trunk@31689 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/doc/news/changes.h | 11 ++- .../include/deal.II/numerics/vector_tools.h | 41 ++++++---- .../deal.II/numerics/vector_tools.templates.h | 78 ++++++++++++++++++- 3 files changed, 111 insertions(+), 19 deletions(-) diff --git a/deal.II/doc/news/changes.h b/deal.II/doc/news/changes.h index 38906d72b7..04671e8690 100644 --- a/deal.II/doc/news/changes.h +++ b/deal.II/doc/news/changes.h @@ -94,6 +94,15 @@ inconvenience this causes.
    +
  1. + 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. +
    + (Wolfgang Bangerth, 2013/11/17) +
    2. +
  2. 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 extract_subvector_to() that allow extracting not just a single value but a whole set.
    - (Fahad Alrasched, 2013/09/02) + (Fahad Alrashed, 2013/09/02)
  3. diff --git a/deal.II/include/deal.II/numerics/vector_tools.h b/deal.II/include/deal.II/numerics/vector_tools.h index d8c6be07ff..47e9823dac 100644 --- a/deal.II/include/deal.II/numerics/vector_tools.h +++ b/deal.II/include/deal.II/numerics/vector_tools.h @@ -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. * * *

    Results

    @@ -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. */ diff --git a/deal.II/include/deal.II/numerics/vector_tools.templates.h b/deal.II/include/deal.II/numerics/vector_tools.templates.h index 7758845376..93e6847b74 100644 --- a/deal.II/include/deal.II/numerics/vector_tools.templates.h +++ b/deal.II/include/deal.II/numerics/vector_tools.templates.h @@ -2404,6 +2404,16 @@ namespace VectorTools }; + template + std::ostream & operator << (std::ostream &out, + const VectorDoFTuple &vdt) + { + for (unsigned int d=0; dcenter() << 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 >::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]; { -- 2.39.5