const unsigned int dofs_per_face = selected_fe->dofs_per_face;
std::vector<unsigned int> dofs_on_face(dofs_per_face);
+
+ // loop over all faces to check
+ // whether they are at a
+ // boundary. note that we need not
+ // take special care of single
+ // lines (using
+ // @p{cell->has_boundary_lines}),
+ // since we do not support
+ // boundaries of dimension dim-2,
+ // and so every boundary line is
+ // also part of a boundary face.
active_face_iterator face = begin_active_face (),
endf = end_face();
for (; face!=endf; ++face)
const unsigned int dofs_per_face = dof.get_fe().dofs_per_face;
std::vector<unsigned int> dofs_on_this_face(dofs_per_face);
- DoFHandler<dim>::active_face_iterator face = dof.begin_active_face(),
- endf = dof.end_face();
+
+ // loop over all faces to check
+ // whether they are at a
+ // boundary. note that we need not
+ // take special care of single
+ // lines (using
+ // @p{cell->has_boundary_lines}),
+ // since we do not support
+ // boundaries of dimension dim-2,
+ // and so every boundary line is
+ // also part of a boundary face.
+ typename DoFHandler<dim>::active_face_iterator face = dof.begin_active_face(),
+ endf = dof.end_face();
for (; face!=endf; ++face)
if (face->at_boundary())
{
DoFHandler<dim>::active_cell_iterator cell = dof.begin_active(),
endc = dof.end();
- (const_cast<Triangulation<dim>& > (dof.get_tria())).clear_user_flags();
+ // clear user flags for further use
+ const_cast<Triangulation<dim>&>(dof.get_tria()).clear_user_flags();
for (; cell!=endc; ++cell)
{
selected_dofs.clear ();
selected_dofs.resize (dof_handler.n_dofs(), false);
std::vector<unsigned int> face_dof_indices (dof_handler.get_fe().dofs_per_face);
+
+ // now loop over all cells and
+ // check whether their faces are at
+ // the boundary. note that we need
+ // not take special care of single
+ // lines being at the boundary
+ // (using
+ // @p{cell->has_boundary_lines}),
+ // since we do not support
+ // boundaries of dimension dim-2,
+ // and so every isolated boundary
+ // line is also part of a boundary
+ // face which we will be visiting
+ // sooner or later
for (DoFHandler<dim>::active_cell_iterator cell=dof_handler.begin_active();
cell!=dof_handler.end(); ++cell)
for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
std::vector<unsigned int> dofs_on_face(dofs_per_face);
unsigned int next_boundary_index = 0;
+ // now loop over all cells and
+ // check whether their faces are at
+ // the boundary. note that we need
+ // not take special care of single
+ // lines being at the boundary
+ // (using
+ // @p{cell->has_boundary_lines}),
+ // since we do not support
+ // boundaries of dimension dim-2,
+ // and so every isolated boundary
+ // line is also part of a boundary
+ // face which we will be visiting
+ // sooner or later
typename DoFHandler<dim>::active_face_iterator face = dof_handler.begin_active_face(),
endf = dof_handler.end_face();
for (; face!=endf; ++face)
template <int dim>
-const StraightBoundary<dim>& Triangulation<dim>::straight_boundary
+const StraightBoundary<dim> & Triangulation<dim>::straight_boundary
= StraightBoundary<dim>();
// interior
for (line_iterator line=begin_line(); line!=end_line(); ++line)
line->set_boundary_indicator (255);
- // next reset all lines bounding boundary
- // quads as on the boundary also
+ // next reset all lines bounding
+ // boundary quads as on the
+ // boundary also. note that since
+ // we are in 3d, there are cases
+ // where one or more lines of a
+ // quad that is not on the
+ // boundary, are actually boundary
+ // lines. they will not be marked
+ // when visiting this
+ // face. however, since we do not
+ // support dim-2 dimensional
+ // boundaries (i.e. internal lines
+ // constituting boundaries), every
+ // such line is also part of a face
+ // that is actually on the
+ // boundary, so sooner or later we
+ // get to mark that line for being
+ // on the boundary
for (quad_iterator quad=begin_quad(); quad!=end_quad(); ++quad)
if (quad->at_boundary())
for (unsigned int l=0; l<4; ++l)
for (cell_iterator cell=begin(0); cell!=end(0); ++cell)
almost_infinite_length += cell->diameter();
- std::vector<double> minimal_length (vertices.size(),
- almost_infinite_length);
+ std::vector<double> minimal_length (vertices.size(),
+ almost_infinite_length);
// also note if a vertex is at
// the boundary
- std::vector<bool> at_boundary (vertices.size(), false);
+ std::vector<bool> at_boundary (vertices.size(), false);
for (active_line_iterator line=begin_active_line();
line != end_line(); ++line)
template <int dim>
void Triangulation<dim>::distort_random (const double factor,
- const bool keep_boundary) {
+ const bool keep_boundary)
+{
// this function is mostly equivalent to
// that for the general dimensional case
// the only difference being the correction
for (cell_iterator cell=begin(0); cell!=end(0); ++cell)
almost_infinite_length += cell->diameter();
- std::vector<double> minimal_length (vertices.size(),
- almost_infinite_length);
+ std::vector<double> minimal_length (vertices.size(),
+ almost_infinite_length);
// also note if a vertex is at
// the boundary
- std::vector<bool> at_boundary (vertices.size(), false);
+ std::vector<bool> at_boundary (vertices.size(), false);
for (active_line_iterator line=begin_active_line();
line != end_line(); ++line)
};
+
template <int dim>
-void Triangulation<dim>::set_all_refine_flags () {
+void Triangulation<dim>::set_all_refine_flags ()
+{
active_cell_iterator cell = begin_active(),
endc = end();
};
+
template <int dim>
-void Triangulation<dim>::refine_global (const unsigned int times) {
+void Triangulation<dim>::refine_global (const unsigned int times)
+{
for (unsigned int i=0; i<times; ++i)
{
set_all_refine_flags();
};
+
// if necessary try to work around a bug in the IBM xlC compiler
#ifdef XLC_WORK_AROUND_STD_BUG
using namespace std;
#if deal_II_dimension == 3
template <>
-void Triangulation<3>::execute_refinement () {
+void Triangulation<3>::execute_refinement ()
+{
const unsigned int dim = 3;
// check whether a new level is needed
vertices[next_unused_vertex]
= boundary[quad->boundary_indicator()]->get_new_point_on_quad (quad);
else
+ // it might be that the
+ // quad itself is not
+ // at the boundary, but
+ // that one of its line
+ // actually is. in this
+ // case, the newly
+ // created vertices at
+ // the centers of the
+ // lines are not
+ // necessarily the mean
+ // values of the
+ // adjacent vertices,
+ // so do not compute
+ // the new vertex as
+ // the mean value of
+ // the 4 vertices of
+ // the face, but rather
+ // as the mean value of
+ // the 8 vertices which
+ // we already have (the
+ // four old ones, and
+ // the four ones
+ // inserted as middle
+ // points for the four
+ // lines). summing up
+ // some more points is
+ // generally cheaper
+ // than first asking
+ // whether one of the
+ // lines is at the
+ // boundary
vertices[next_unused_vertex]
= (quad->vertex(0) + quad->vertex(1) +
- quad->vertex(2) + quad->vertex(3)) / 4;
+ quad->vertex(2) + quad->vertex(3) +
+ quad->line(0)->child(0)->vertex(1) +
+ quad->line(1)->child(0)->vertex(1) +
+ quad->line(2)->child(0)->vertex(1) +
+ quad->line(3)->child(0)->vertex(1) ) / 8;
// now that we created the right
// point, make up the four
std::map<unsigned int,double> boundary_values;
if (enforce_zero_boundary == true)
- // no need to project boundary values, but
- // enforce homogeneous boundary values
+ // no need to project boundary
+ // values, but enforce
+ // homogeneous boundary values
// anyway
{
- DoFHandler<dim>::active_face_iterator face = dof.begin_active_face(),
- endf = dof.end_face();
+ // loop over all boundary faces
+ // to get all dof indices of
+ // dofs on the boundary. note
+ // that in 3d there are cases
+ // where a face is not at the
+ // boundary, yet one of its
+ // lines is, and we should
+ // consider the degrees of
+ // freedom on it as boundary
+ // nodes. likewise, in 2d and
+ // 3d there are cases where a
+ // cell is only at the boundary
+ // by one vertex. nevertheless,
+ // since we do not support
+ // boundaries with dimension
+ // less or equal to dim-2, each
+ // such boundary dof is also
+ // found from some other face
+ // that is actually wholly on
+ // the boundary, not only by
+ // one line or one vertex
+ typename DoFHandler<dim>::active_face_iterator face = dof.begin_active_face(),
+ endf = dof.end_face();
std::vector<unsigned int> face_dof_indices (fe.dofs_per_face);
for (; face!=endf; ++face)
if (face->at_boundary())