<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
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>
/**
- * 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
* 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>
* 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.
*/
* 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.
*/
* 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.
*/
* 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.
*/
* "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.
*/
/**
* 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.
*/
* 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.
*/
* 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.
*/
* 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.
*/
* 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.
*/
};
+ 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
.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
}
}
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
= 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
// 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:
{
// 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
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];
{
projects / dealii.git / commitdiff