unsigned int number_of_children () const;
/**
- * Enquire some information about
- * whether the face with the
- * given number is in standard
- * orientation or not. This
- * information is only useful
- * (and available) in 3d, see the
+ * Return whether the face with
+ * index @p{face} has its normal
+ * pointing in the standard
+ * direction (@p{true}) or
+ * whether it is the opposite
+ * (@p{false}). Which is the
+ * standard direction is
+ * documented with the
+ * @ref{Triangulation} class. In
+ * 1d and 2d, this is always
+ * @p{true}, but in 3d it may be
+ * different, see the respective
+ * discussion in the
* documentation of the
- * respective class. For all
- * other dimensions, this
- * function should not be called
- * and triggers and assertion if
- * done so.
+ * @ref{Triangulation} class.
+ *
+ * This function is really only
+ * for internal use in the
+ * library unless you absolutely
+ * know what this is all about.
*/
- bool get_face_orientation (const unsigned int face) const;
+ bool face_orientation (const unsigned int face) const;
private:
/**
unsigned int number_of_children () const;
/**
- * Enquire some information about
- * whether the face with the
- * given number is in standard
- * orientation or not. This
- * information is only useful
- * (and available) in 3d, see the
+ * Return whether the face with
+ * index @p{face} has its normal
+ * pointing in the standard
+ * direction (@p{true}) or
+ * whether it is the opposite
+ * (@p{false}). Which is the
+ * standard direction is
+ * documented with the
+ * @ref{Triangulation} class. In
+ * 1d and 2d, this is always
+ * @p{true}, but in 3d it may be
+ * different, see the respective
+ * discussion in the
* documentation of the
- * respective class. For all
- * other dimensions, this
- * function should not be called
- * and triggers and assertion if
- * done so.
+ * @ref{Triangulation} class.
+ *
+ * This function is really only
+ * for internal use in the
+ * library unless you absolutely
+ * know what this is all about.
*/
- bool get_face_orientation (const unsigned int face) const;
+ bool face_orientation (const unsigned int face) const;
private:
/**
unsigned int number_of_children () const;
/**
- * Enquire some information about
- * whether the face with the
- * given number is in standard
- * orientation or not. This
- * information is only useful
- * (and available) in 3d, see the
+ * Return whether the face with
+ * index @p{face} has its normal
+ * pointing in the standard
+ * direction (@p{true}) or
+ * whether it is the opposite
+ * (@p{false}). Which is the
+ * standard direction is
+ * documented with the
+ * @ref{Triangulation} class. In
+ * 1d and 2d, this is always
+ * @p{true}, but in 3d it may be
+ * different, see the respective
+ * discussion in the
* documentation of the
- * respective class. For all
- * other dimensions, this
- * function should not be called
- * and triggers and assertion if
- * done so.
+ * @ref{Triangulation} class.
+ *
+ * This function is really only
+ * for internal use in the
+ * library unless you absolutely
+ * know what this is all about.
*/
- bool get_face_orientation (const unsigned int face) const;
+ bool face_orientation (const unsigned int face) const;
private:
/**
unsigned int number_of_children () const;
/**
- * Return whether the quad with
+ * Return whether the face with
* index @p{face} has its normal
* pointing in the standard
* direction (@p{true}) or
* (@p{false}). Which is the
* standard direction is
* documented with the
+ * @ref{Triangulation} class. In
+ * 1d and 2d, this is always
+ * @p{true}, but in 3d it may be
+ * different, see the respective
+ * discussion in the
+ * documentation of the
* @ref{Triangulation} class.
*
* This function is really only
* library unless you absolutely
* know what this is all about.
*/
- bool get_face_orientation (const unsigned int face) const;
+ bool face_orientation (const unsigned int face) const;
/**
* Set whether the quad with
#include <cmath>
+//TODO[WB]: TriaObjectAccessor<dim,N> should never be used since we have specializations. Can we remove the implementations of functions, or what are they there for?
+
/*------------------------ Functions: TriaAccessor ---------------------------*/
+template <int dim>
+inline
+bool
+TriaObjectAccessor<1,dim>::face_orientation (const unsigned int) const
+{
+ return true;
+}
+
+
+
template <int dim>
inline
void
+template <int dim>
+inline
+bool
+TriaObjectAccessor<2,dim>::face_orientation (const unsigned int) const
+{
+ return true;
+}
+
+
+
template <int dim>
inline
void
{ 4, 3, 0 }};
return (this->quad(lookup_table[i][0])
- ->line(get_face_orientation(lookup_table[i][0]) ?
+ ->line(face_orientation(lookup_table[i][0]) ?
lookup_table[i][1] :
lookup_table[i][2]));
}
{ 4, 3, 0 }};
return (this->quad(lookup_table[i][0])
- ->line_index(get_face_orientation(lookup_table[i][0]) ?
+ ->line_index(face_orientation(lookup_table[i][0]) ?
lookup_table[i][1] :
lookup_table[i][2]));
}
+template <int dim>
+inline
+bool
+TriaObjectAccessor<3, dim>::
+face_orientation (const unsigned int face) const
+{
+ Assert (used(), typename TriaAccessor<dim>::ExcCellNotUsed());
+ Assert (face<GeometryInfo<3>::faces_per_cell,
+ ExcIndexRange (face, 0, GeometryInfo<3>::faces_per_cell));
+ Assert (this->present_index * GeometryInfo<3>::faces_per_cell + face
+ < this->tria->levels[this->present_level]
+ ->hexes.face_orientations.size(),
+ ExcInternalError());
+
+ return (this->tria->levels[this->present_level]
+ ->hexes.face_orientations[this->present_index *
+ GeometryInfo<3>::faces_per_cell
+ +
+ face]);
+}
+
+
+
template <int dim>
inline
void
// face or subface
case 3:
return ((face_no +
- (cell->get_face_orientation(face_no) == true ?
+ (cell->face_orientation(face_no) == true ?
0 : GeometryInfo<dim>::faces_per_cell))
* n_quadrature_points);
Assert (subface_no+1 < GeometryInfo<dim>::subfaces_per_face+1,
ExcInternalError());
Assert (n_quadrature_points > 0, ExcInternalError());
+
+ Assert (cell->has_children() == false,
+ ExcMessage ("You can't use subface data for cells that are "
+ "already refined. Iterate over their children "
+ "instead in these cases."));
switch (dim)
{
case 3:
return (((face_no * GeometryInfo<dim>::subfaces_per_face +
subface_no)
- + (cell->get_face_orientation(face_no) == true ?
+ + (cell->face_orientation(face_no) == true ?
0 :
GeometryInfo<dim>::faces_per_cell *
- GeometryInfo<dim>::subfaces_per_face))
+ GeometryInfo<dim>::subfaces_per_face)
+ )
* n_quadrature_points);
default:
Assert (false, ExcInternalError());
// opposite direction, so
// store that up front
const bool face_orientation[6]
- = { hex->get_face_orientation (0),
- hex->get_face_orientation (1),
- hex->get_face_orientation (2),
- hex->get_face_orientation (3),
- hex->get_face_orientation (4),
- hex->get_face_orientation (5) };
+ = { hex->face_orientation (0),
+ hex->face_orientation (1),
+ hex->face_orientation (2),
+ hex->face_orientation (3),
+ hex->face_orientation (4),
+ hex->face_orientation (5) };
const unsigned int line_indices[30]
= {
for (unsigned int f=0; f<GeometryInfo<dim>::faces_per_cell; ++f)
for (unsigned int s=0; s<GeometryInfo<dim>::subfaces_per_face; ++s)
new_hexes[GeometryInfo<dim>::child_cell_on_face(f,s)]
- ->set_face_orientation(f, hex->get_face_orientation(f));
+ ->set_face_orientation(f, hex->face_orientation(f));
/////////////////////////////////
// all
// wrong...)
const bool orient
- = (neighbor->get_face_orientation(nb_nb)
+ = (neighbor->face_orientation(nb_nb)
&&
- hex->get_face_orientation(face));
+ hex->face_orientation(face));
static const unsigned int
child_switch_table[GeometryInfo<dim>::subfaces_per_face]
-template <int dim>
-bool
-TriaObjectAccessor<1,dim>::get_face_orientation (const unsigned int) const
-{
- Assert (false, ExcInternalError());
- return true;
-}
-
-
/*------------------------ Functions: QuadAccessor ---------------------------*/
template <int dim>
-template <int dim>
-bool
-TriaObjectAccessor<2,dim>::get_face_orientation (const unsigned int) const
-{
- Assert (false, ExcInternalError());
- return true;
-}
-
-
/*------------------------ Functions: TriaObjectAccessor ---------------------------*/
{ 0, 3, 2, 1 };
if (corner<4)
{
- if (get_face_orientation(0) == true)
+ if (face_orientation(0) == true)
return quad(0)->vertex_index(corner);
else
return quad(0)->vertex_index(vertex_translation[corner]);
}
else
{
- if (get_face_orientation(1) == true)
+ if (face_orientation(1) == true)
return quad(1)->vertex_index(corner-4);
else
return quad(1)->vertex_index(vertex_translation[corner-4]);
-template <int dim>
-bool
-TriaObjectAccessor<3, dim>::
-get_face_orientation (const unsigned int face) const
-{
- Assert (used(), typename TriaAccessor<dim>::ExcCellNotUsed());
- Assert (face<GeometryInfo<3>::faces_per_cell,
- ExcIndexRange (face, 0, GeometryInfo<3>::faces_per_cell));
- Assert (this->present_index * GeometryInfo<3>::faces_per_cell + face
- < this->tria->levels[this->present_level]
- ->hexes.face_orientations.size(),
- ExcInternalError());
-
- return (this->tria->levels[this->present_level]
- ->hexes.face_orientations[this->present_index *
- GeometryInfo<3>::faces_per_cell
- +
- face]);
-}
-
-
-
template <int dim>
void
TriaObjectAccessor<3, dim>::
{
const unsigned int n_solution_vectors = solutions.size();
- // make up a fe face values object for the
- // restriction of the finite element function
- // to a face, for the present cell and its
- // neighbor. In principle we would only need
- // one at a time, but this way we can
- // have more fine grained access to what
- // values really need to be computed (we
- // need not compute all values on the
- // neighbor cells, so using two objects
+ // make up a fe face values object
+ // for the restriction of the
+ // finite element function to a
+ // face, for the present cell and
+ // its neighbor. In principle we
+ // would only need one at a time,
+ // but this way we can have more
+ // fine grained access to what
+ // values really need to be
+ // computed (we need not compute
+ // all values on the neighbor
+ // cells, so using two objects
// gives us a performance gain).
+ //
+ // in debug mode, make sure that
+ // some data matches, so compute
+ // quadrature points always
FEFaceValues<dim> fe_face_values_cell (mapping,
dof_handler.get_fe(),
quadrature,
- UpdateFlags(update_gradients |
- update_JxW_values |
- ((!neumann_bc.empty() ||
- (coefficients != 0)) ?
- update_q_points : 0) |
- update_normal_vectors));
+ UpdateFlags(
+ update_gradients |
+ update_JxW_values |
+ ((!neumann_bc.empty() ||
+ (coefficients != 0)) ?
+ update_q_points : 0) |
+#ifdef DEBUG
+ update_q_points |
+#endif
+ update_normal_vectors));
FEFaceValues<dim> fe_face_values_neighbor (mapping,
dof_handler.get_fe(),
quadrature,
FESubfaceValues<dim> fe_subface_values (mapping,
dof_handler.get_fe(),
quadrature,
- update_gradients);
+ UpdateFlags(
+ update_gradients
+#ifdef DEBUG
+ | update_q_points
+#endif
+ ));
active_cell_iterator cell = dof_handler.begin_active();
{
// get an iterator pointing to the
// cell behind the present subface
- const bool orientation_flag
- = (dim != 3 ?
- true :
- cell->get_face_orientation(face_no) == true
- &&
- neighbor->get_face_orientation(neighbor_neighbor) == true);
static const unsigned int subface_translation[4]
= { 0, 3, 2, 1 };
+ // see whether face and
+ // the neighbor's
+ // counterface share the
+ // same indexing of
+ // children. if not so,
+ // translate child
+ // indices
+ const bool face_orientations_match
+ = (neighbor->face_orientation(neighbor_neighbor) ==
+ cell->face_orientation(face_no));
const unsigned int neighbor_child_index
= (GeometryInfo<dim>::
child_cell_on_face(neighbor_neighbor,
- (orientation_flag ?
+ (face_orientations_match ?
subface_no :
subface_translation[subface_no])));
const active_cell_iterator neighbor_child
= neighbor->child(neighbor_child_index);
- Assert (neighbor_child->face(neighbor_neighbor) ==
- cell->face(face_no)->child(subface_no),
- ExcInternalError());
Assert (!neighbor->child(neighbor_child_index)->has_children(),
ExcInternalError());
- // restrict the finite element on the
- // present cell to the subface and
- // store the gradient of the solution
- // in psi
+ // restrict the finite element
+ // on the present cell to the
+ // subface
fe_subface_values.reinit (cell, face_no, subface_no);
+ // restrict the finite element
+ // on the neighbor cell to the
+ // common @p{subface}.
+ fe_face_values.reinit (neighbor_child, neighbor_neighbor);
+
+ // make sure that quadrature
+ // points match
+ for (unsigned int q=0; q<n_q_points; ++q)
+ Assert ((fe_face_values.quadrature_point(q) -
+ fe_subface_values.quadrature_point(q)).square()
+ <
+ 1.e-15 * (fe_face_values.quadrature_point(q).square() +
+ fe_subface_values.quadrature_point(q).square()),
+ ExcInternalError());
+
+ // store the gradient of the
+ // solution in psi
for (unsigned int n=0; n<n_solution_vectors; ++n)
fe_subface_values.get_function_grads (*solutions[n], per_thread_data.psi[n]);
- // restrict the finite element on the
- // neighbor cell to the common @p{subface}.
- // store the gradient in @p{neighbor_psi}
-
- fe_face_values.reinit (neighbor_child, neighbor_neighbor);
-
+ // store the gradient from the
+ // neighbor's side in
+ // @p{neighbor_psi}
for (unsigned int n=0; n<n_solution_vectors; ++n)
fe_face_values.get_function_grads (*solutions[n], per_thread_data.neighbor_psi[n]);
if (coefficients->n_components == 1)
{
coefficients->value_list (fe_face_values.get_quadrature_points(),
- per_thread_data.coefficient_values1);
+ per_thread_data.coefficient_values1);
for (unsigned int n=0; n<n_solution_vectors; ++n)
for (unsigned int component=0; component<n_components; ++component)
for (unsigned int point=0; point<n_q_points; ++point)