/**
* This class is a helper class to facilitate the usage of quadrature formulae
* on faces or subfaces of cells. It computes the locations of quadrature
- * points on the unit cell from a quadrature object for a mannifold of
+ * points on the unit cell from a quadrature object for a manifold of
* one dimension less than that of the cell and the number of the face.
* For example, giving the Simpson rule in one dimension and using the
* @p{project_to_face} function with face number 1, the returned points will
* unit cell), the same applies as above. Note the order in which the
* children of a face are numbered, which in two dimensions coincides
* with the orientation of the face.
+ *
+ * The second set of functions generates a quadrature formula by
+ * projecting a given quadrature rule on @em{all} faces and
+ * subfaces. This is used in the @ref{FEFaceValues} and
+ * @ref{FESubfaceValues} classes. Since we now have the quadrature
+ * points of all faces and subfaces in one array, we need to have a
+ * way to find the starting index of the points and weights
+ * corresponding to one face or subface within this array. This is
+ * done through the @ref{DataSetDescriptor} member class.
*
- * The different functions are grouped into a common class to avoid putting
- * them into global namespace (and to make documentation easier, since
- * presently the documentation tool can only handle classes, not global
- * functions). However, since they have no local data, all functions are
- * declared @p{static} and can be called without creating an object of this
- * class.
+ * The different functions are grouped into a common class to avoid
+ * putting them into global namespace. However, since they have no
+ * local data, all functions are declared @p{static} and can be
+ * called without creating an object of this class.
+ *
+ * For the 3d case, you should note that the orientation of faces is
+ * even more intricate than for two dimensions. Quadrature formulae
+ * are projected upon the faces in their standard orientation, not to
+ * the inside or outside of the hexahedron. Refer to the
+ * documentation of the @p{Triangulation} class for a description of
+ * the orientation of the different faces. To make things more
+ * complicated, in 3d we allow faces in two orientations (which can
+ * be identified using @p{cell->face_orientation(face)}), so we have
+ * to project quadrature formula onto faces and subfaces in two
+ * orientations. The @ref{DataSetDescriptor} member class is used to
+ * identify where each dataset starts.
*
- * For the 3d case, you should note that the orientation of faces is even
- * more intricate than for two dimensions. Quadrature formulae are projected
- * upon the faces in their standard orientation, not to the inside or outside
- * of the hexahedron! Refer to the documentation of the @p{Triangulation} class
- * for a description of the orientation of the different faces.
+ * @author Wolfgang Bangerth, 1998, 1999, 2003
*/
template <int dim>
class QProjector
project_to_child (const Quadrature<dim> &quadrature,
const unsigned int child_no);
+ /**
+ * Since the
+ * @p{project_to_all_faces} and
+ * @p{project_to_all_subfaces}
+ * functions chain together the
+ * quadrature points and weights
+ * of all projections of a face
+ * quadrature formula to the
+ * faces or subfaces of a cell,
+ * we need a way to identify
+ * where the starting index of
+ * the points and weights for a
+ * particular face or subface
+ * is. This class provides this:
+ * there are static member
+ * functions that generate
+ * objects of this type, given
+ * face or subface indices, and
+ * you can then use the generated
+ * object in place of an integer
+ * that denotes the offset of a
+ * given dataset.
+ *
+ * @author Wolfgang Bangerth, 2003
+ */
+ class DataSetDescriptor
+ {
+ public:
+ /**
+ * Default constructor. This
+ * doesn't do much except
+ * generating an invalid
+ * index, since you didn't
+ * give a valid descriptor of
+ * the cell, face, or subface
+ * you wanted.
+ */
+ DataSetDescriptor ();
+
+ /**
+ * Static function to
+ * generate the offset of a
+ * cell. Since we only have
+ * one cell per quadrature
+ * object, this offset is of
+ * course zero, but we carry
+ * this function around for
+ * consistency with the other
+ * static functions.
+ */
+ static DataSetDescriptor cell ();
+
+ /**
+ * Static function to
+ * generate an offset object
+ * for a given face of a cell
+ * with the given face
+ * orientation. This function
+ * of course is only allowed
+ * if @p{dim>=2}, and the
+ * face orientation is
+ * ignored if the space
+ * dimension equals 2.
+ *
+ * The last argument denotes
+ * the number of quadrature
+ * points the
+ * lower-dimensional face
+ * quadrature formula (the
+ * one that has been
+ * projected onto the faces)
+ * has.
+ */
+ static
+ DataSetDescriptor
+ face (const unsigned int face_no,
+ const bool face_orientation,
+ const unsigned int n_quadrature_points);
+
+ /**
+ * Static function to
+ * generate an offset object
+ * for a given subface of a
+ * cell with the given face
+ * orientation. This function
+ * of course is only allowed
+ * if @p{dim>=2}, and the
+ * face orientation is
+ * ignored if the space
+ * dimension equals 2.
+ *
+ * The last argument denotes
+ * the number of quadrature
+ * points the
+ * lower-dimensional face
+ * quadrature formula (the
+ * one that has been
+ * projected onto the faces)
+ * has.
+ */
+ static
+ DataSetDescriptor
+ sub_face (const unsigned int face_no,
+ const unsigned int subface_no,
+ const bool face_orientation,
+ const unsigned int n_quadrature_points);
+
+ /**
+ * Conversion operator to an
+ * integer denoting the
+ * offset of the first
+ * element of this dataset in
+ * the set of quadrature
+ * formulas all projected
+ * onto faces and
+ * subfaces. This conversion
+ * operator allows us to use
+ * offset descriptor objects
+ * in place of integer
+ * offsets.
+ */
+ operator unsigned int () const;
+
+ private:
+ /**
+ * Store the integer offset
+ * for a given cell, face, or
+ * subface.
+ */
+ const unsigned int dataset_offset;
+
+ /**
+ * This is the real
+ * constructor, but it is
+ * private and thus only
+ * available to the static
+ * member functions above.
+ */
+ DataSetDescriptor (const unsigned int dataset_offset);
+ };
+
private:
/**
* Given a quadrature object in
+template <int dim>
+typename QProjector<dim>::DataSetDescriptor
+QProjector<dim>::DataSetDescriptor::cell ()
+{
+ return 0;
+}
+
+
+template <int dim>
+typename QProjector<dim>::DataSetDescriptor
+QProjector<dim>::DataSetDescriptor::
+face (const unsigned int face_no,
+ const bool face_orientation,
+ const unsigned int n_quadrature_points)
+{
+ Assert (dim != 1, ExcInternalError());
+ Assert (face_no < GeometryInfo<dim>::faces_per_cell,
+ ExcInternalError());
+ Assert (n_quadrature_points > 0, ExcInternalError());
+
+ switch (dim)
+ {
+ case 1:
+ case 2:
+ return face_no * n_quadrature_points;
+
+ // in 3d, we have to
+ // account for faces
+ // that have reverse
+ // orientation. thus,
+ // we have to store
+ // _two_ data sets per
+ // face or subface
+ case 3:
+ return ((face_no +
+ (face_orientation == true ?
+ 0 : GeometryInfo<dim>::faces_per_cell))
+ * n_quadrature_points);
+
+ default:
+ Assert (false, ExcInternalError());
+ }
+ return static_cast<unsigned int>(-1);
+}
+
+
+
+template <int dim>
+typename QProjector<dim>::DataSetDescriptor
+QProjector<dim>::DataSetDescriptor::
+sub_face (const unsigned int face_no,
+ const unsigned int subface_no,
+ const bool face_orientation,
+ const unsigned int n_quadrature_points)
+{
+ Assert (dim != 1, ExcInternalError());
+ Assert (face_no < GeometryInfo<dim>::faces_per_cell,
+ ExcInternalError());
+ // the trick with +1 prevents
+ // that we get a warning in 1d
+ Assert (subface_no+1 < GeometryInfo<dim>::subfaces_per_face+1,
+ ExcInternalError());
+ Assert (n_quadrature_points > 0, ExcInternalError());
+
+ switch (dim)
+ {
+ case 1:
+ case 2:
+ return ((face_no * GeometryInfo<dim>::subfaces_per_face +
+ subface_no)
+ * n_quadrature_points);
+
+ // for 3d, same as above:
+ case 3:
+ return (((face_no * GeometryInfo<dim>::subfaces_per_face +
+ subface_no)
+ + (face_orientation == true ?
+ 0 :
+ GeometryInfo<dim>::faces_per_cell *
+ GeometryInfo<dim>::subfaces_per_face)
+ )
+ * n_quadrature_points);
+ default:
+ Assert (false, ExcInternalError());
+ }
+ return static_cast<unsigned int>(-1);
+}
+
+
+template <int dim>
+QProjector<dim>::DataSetDescriptor::operator unsigned int () const
+{
+ return dataset_offset;
+}
+
+
+
+template <int dim>
+QProjector<dim>::DataSetDescriptor::
+DataSetDescriptor (const unsigned int dataset_offset)
+ :
+ dataset_offset (dataset_offset)
+{}
+
+
+template <int dim>
+QProjector<dim>::DataSetDescriptor::
+DataSetDescriptor ()
+ :
+ dataset_offset (static_cast<unsigned int>(-1))
+{}
+
+
+
// ------------------------------------------------------------ //
#include <base/config.h>
#include <base/table.h>
-#include <cmath>
+#include <base/quadrature.h>
#include <grid/tria_iterator.h>
#include <dofs/dof_accessor.h>
#include <fe/mapping.h>
+#include <cmath>
-/*!@addtogroup fe */
-/*@{*/
-namespace internal
-{
- template <int dim>
- class DataSetDescriptor
- {
- public:
- DataSetDescriptor ();
-
- static DataSetDescriptor cell ();
-
- static
- DataSetDescriptor
- face (const typename DoFHandler<dim>::cell_iterator &cell,
- const unsigned int face_no,
- const unsigned int n_quadrature_points);
-
- static
- DataSetDescriptor
- sub_face (const typename DoFHandler<dim>::cell_iterator &cell,
- const unsigned int face_no,
- const unsigned int subface_no,
- const unsigned int n_quadrature_points);
-
- unsigned int offset () const;
-
- private:
- const unsigned int dataset_offset;
-
- DataSetDescriptor (const unsigned int dataset_offset);
- };
-}
+/*!@addtogroup fe */
+/*@{*/
+
+
/**
* Mapping of general quadrilateral/hexahedra by d-linear shape
* functions.
DeclException0 (ExcAccessToUninitializedField);
protected:
+
/**
- * Typedef the right data set
- * descriptor to a local type to
- * make use somewhat simpler.
+ * Declare a convenience typedef
+ * for the class that describes
+ * offsets into quadrature
+ * formulas projected onto faces
+ * and subfaces.
*/
- typedef internal::DataSetDescriptor<dim> DataSetDescriptor;
+ typedef
+ typename QProjector<dim>::DataSetDescriptor
+ DataSetDescriptor;
/**
* Implementation of the interface in
* 3 and 2, etc.
*
* For three spatial dimensions,
- * the exact order of the children is
- * laid down in the documentation of
- * the @ref{Triangulation} class.
+ * the exact order of the
+ * children is laid down in the
+ * documentation of the
+ * @ref{Triangulation} class.
+ * However, it must be noted that
+ * this class and function only
+ * deals with faces in standard
+ * orientation. In 3d, faces can
+ * exist in two orientations,
+ * though, and if a face is in
+ * the wrong orientation, then
+ * this function may not give you
+ * what you want. You can inquire
+ * about the face orientation
+ * using the
+ * @p{cell->face_orientation}
+ * function, and the function to
+ * ask for a neighbor's cell
+ * behind a given face and
+ * subface is
+ * @p{cell->neighbor_child_on_subface}.
+ * The latter function, in
+ * contrast to the present one,
+ * also takes into account the
+ * actual orientation of the
+ * faces of a cell and will
+ * return the correct result in
+ * all cases.
*/
static unsigned int child_cell_on_face (const unsigned int face,
const unsigned int subface);
* 3 and 2, etc.
*
* For three spatial dimensions,
- * the exact order of the children is
- * laid down in the documentation of
- * the @ref{Triangulation} class.
+ * the exact order of the
+ * children is laid down in the
+ * documentation of the
+ * @ref{Triangulation}
+ * class. However, it must be
+ * noted that this class and
+ * function only deals with faces
+ * in standard orientation. In
+ * 3d, faces can exist in two
+ * orientations, though, and if a
+ * face is in the wrong
+ * orientation, then this
+ * function may not give you what
+ * you want. You can inquire
+ * about the face orientation
+ * using the
+ * @p{cell->face_orientation}
+ * function, and the function to
+ * ask for a neighbor's cell
+ * behind a given face and
+ * subface is
+ * @p{cell->neighbor_child_on_subface}.
+ * The latter function, in
+ * contrast to the present one,
+ * also takes into account the
+ * actual orientation of the
+ * faces of a cell and will
+ * return the correct result in
+ * all cases.
*/
static unsigned int child_cell_on_face (const unsigned int face,
const unsigned int subface);
* 3 and 2, etc.
*
* For three spatial dimensions,
- * the exact order of the children is
- * laid down in the documentation of
- * the @ref{Triangulation} class.
+ * the exact order of the
+ * children is laid down in the
+ * documentation of the
+ * @ref{Triangulation}
+ * class. However, it must be
+ * noted that this class and
+ * function only deals with faces
+ * in standard orientation. In
+ * 3d, faces can exist in two
+ * orientations, though, and if a
+ * face is in the wrong
+ * orientation, then this
+ * function may not give you what
+ * you want. You can inquire
+ * about the face orientation
+ * using the
+ * @p{cell->face_orientation}
+ * function, and the function to
+ * ask for a neighbor's cell
+ * behind a given face and
+ * subface is
+ * @p{cell->neighbor_child_on_subface}.
+ * The latter function, in
+ * contrast to the present one,
+ * also takes into account the
+ * actual orientation of the
+ * faces of a cell and will
+ * return the correct result in
+ * all cases.
*/
static unsigned int child_cell_on_face (const unsigned int face,
const unsigned int subface);
* 3 and 2, etc.
*
* For three spatial dimensions,
- * the exact order of the children is
- * laid down in the documentation of
- * the @ref{Triangulation} class.
+ * the exact order of the
+ * children is laid down in the
+ * documentation of the
+ * @ref{Triangulation}
+ * class. However, it must be
+ * noted that this class and
+ * function only deals with faces
+ * in standard orientation. In
+ * 3d, faces can exist in two
+ * orientations, though, and if a
+ * face is in the wrong
+ * orientation, then this
+ * function may not give you what
+ * you want. You can inquire
+ * about the face orientation
+ * using the
+ * @p{cell->face_orientation}
+ * function, and the function to
+ * ask for a neighbor's cell
+ * behind a given face and
+ * subface is
+ * @p{cell->neighbor_child_on_subface}.
+ * The latter function, in contrast
+ * to the present one, also takes
+ * into account the actual
+ * orientation of the faces of a
+ * cell and will return the
+ * correct result in all cases.
*/
static unsigned int child_cell_on_face (const unsigned int face,
const unsigned int subface);
if (neighbor->level() < cell->level())
continue;
- unsigned int neighbor_face = cell->neighbor_of_neighbor(face);
+ const unsigned int neighbor_face
+ = cell->neighbor_of_neighbor(face);
if (neighbor->has_children())
{
sub_nr != GeometryInfo<dim>::subfaces_per_face;
++sub_nr)
{
- typename DoFHandler<dim>::cell_iterator sub_neighbor
- = neighbor->
- child(GeometryInfo<dim>::child_cell_on_face(neighbor_face, sub_nr));
+ const typename DoFHandler<dim>::cell_iterator
+ sub_neighbor
+ = cell->neighbor_child_on_subface (face, sub_nr);
sub_neighbor->get_dof_indices (dofs_on_other_cell);
for (unsigned int i=0; i<total_dofs; ++i)
- {
- for (unsigned int j=0; j<total_dofs; ++j)
- {
- sparsity.add (dofs_on_this_cell[i],
- dofs_on_other_cell[j]);
- sparsity.add (dofs_on_other_cell[i],
- dofs_on_this_cell[j]);
- }
- }
+ for (unsigned int j=0; j<total_dofs; ++j)
+ {
+ sparsity.add (dofs_on_this_cell[i],
+ dofs_on_other_cell[j]);
+ sparsity.add (dofs_on_other_cell[i],
+ dofs_on_this_cell[j]);
+ }
+
sub_neighbor->face(neighbor_face)->set_user_flag ();
}
- } else {
+ }
+ else
+ {
neighbor->get_dof_indices (dofs_on_other_cell);
for (unsigned int i=0; i<total_dofs; ++i)
- {
- for (unsigned int j=0; j<total_dofs; ++j)
- {
- sparsity.add (dofs_on_this_cell[i],
- dofs_on_other_cell[j]);
- sparsity.add (dofs_on_other_cell[i],
- dofs_on_this_cell[j]);
- }
- }
+ for (unsigned int j=0; j<total_dofs; ++j)
+ {
+ sparsity.add (dofs_on_this_cell[i],
+ dofs_on_other_cell[j]);
+ sparsity.add (dofs_on_other_cell[i],
+ dofs_on_this_cell[j]);
+ }
+
neighbor->face(neighbor_face)->set_user_flag ();
}
}
face < GeometryInfo<dim>::faces_per_cell;
++face)
{
- typename DoFHandler<dim>::face_iterator cell_face = cell->face(face);
+ const typename DoFHandler<dim>::face_iterator
+ cell_face = cell->face(face);
if (cell_face->user_flag_set ())
continue;
const bool j_non_zero_i = support_on_face (j, face);
if (flux_dof_mask(i,j) == 'f')
- {
- sparsity.add (dofs_on_this_cell[i],
- dofs_on_this_cell[j]);
- }
+ sparsity.add (dofs_on_this_cell[i],
+ dofs_on_this_cell[j]);
if (flux_dof_mask(i,j) == 'a')
- {
- if (i_non_zero_i && j_non_zero_i)
- sparsity.add (dofs_on_this_cell[i],
- dofs_on_this_cell[j]);
- }
+ if (i_non_zero_i && j_non_zero_i)
+ sparsity.add (dofs_on_this_cell[i],
+ dofs_on_this_cell[j]);
}
}
}
if (neighbor->level() < cell->level())
continue;
- unsigned int neighbor_face = cell->neighbor_of_neighbor(face);
+ const unsigned int
+ neighbor_face = cell->neighbor_of_neighbor(face);
if (neighbor->has_children())
{
sub_nr != GeometryInfo<dim>::subfaces_per_face;
++sub_nr)
{
- typename DoFHandler<dim>::cell_iterator sub_neighbor
- = neighbor->
- child(GeometryInfo<dim>::child_cell_on_face(neighbor_face, sub_nr));
+ const typename DoFHandler<dim>::cell_iterator
+ sub_neighbor
+ = cell->neighbor_child_on_subface (face, sub_nr);
sub_neighbor->get_dof_indices (dofs_on_other_cell);
for (unsigned int i=0; i<total_dofs; ++i)
}
sub_neighbor->face(neighbor_face)->set_user_flag ();
}
- } else {
+ }
+ else
+ {
neighbor->get_dof_indices (dofs_on_other_cell);
for (unsigned int i=0; i<total_dofs; ++i)
{
#include <fe/fe.h>
#include <fe/mapping.h>
#include <fe/fe_dgp.h>
-#include <fe/mapping_q1.h>
#include <fe/fe_values.h>
#ifdef HAVE_STD_STRINGSTREAM
if (flags & update_second_derivatives)
this->compute_2nd (mapping, cell,
- internal::DataSetDescriptor<dim>::cell().offset(),
+ QProjector<dim>::DataSetDescriptor::cell(),
mapping_data, fe_data, data);
}
// offset determines which data set
// to take (all data sets for all
// faces are stored contiguously)
- const unsigned int offset
- = (internal::DataSetDescriptor<dim>::
- face (cell, face,
- quadrature.n_quadrature_points)).offset();
+ const typename QProjector<dim>::DataSetDescriptor offset
+ = (QProjector<dim>::DataSetDescriptor::
+ face (face, cell->face_orientation(face),
+ quadrature.n_quadrature_points));
const UpdateFlags flags(fe_data.update_once | fe_data.update_each);
// offset determines which data set
// to take (all data sets for all
// sub-faces are stored contiguously)
- const unsigned int offset
- = (internal::DataSetDescriptor<dim>::
- sub_face (cell, face, subface,
- quadrature.n_quadrature_points)).offset();
+ const typename QProjector<dim>::DataSetDescriptor offset
+ = (QProjector<dim>::DataSetDescriptor::
+ sub_face (face, subface, cell->face_orientation(face),
+ quadrature.n_quadrature_points));
const UpdateFlags flags(fe_data.update_once | fe_data.update_each);
#include <grid/tria_iterator.h>
#include <dofs/dof_accessor.h>
#include <fe/fe.h>
-#include <fe/mapping_q1.h>
#include <fe/mapping.h>
#include <fe/fe_dgq.h>
#include <fe/fe_values.h>
if (flags & update_second_derivatives)
this->compute_2nd (mapping, cell,
- internal::DataSetDescriptor<dim>::cell().offset(),
+ QProjector<dim>::DataSetDescriptor::cell(),
mapping_data, fe_data, data);
}
// offset determines which data set
// to take (all data sets for all
// faces are stored contiguously)
- const unsigned int offset
- = (internal::DataSetDescriptor<dim>::
- face (cell, face,
- quadrature.n_quadrature_points)).offset();
+ const typename QProjector<dim>::DataSetDescriptor offset
+ = (QProjector<dim>::DataSetDescriptor::
+ face (face, cell->face_orientation(face),
+ quadrature.n_quadrature_points));
const UpdateFlags flags(fe_data.update_once | fe_data.update_each);
// offset determines which data set
// to take (all data sets for all
// sub-faces are stored contiguously)
- const unsigned int offset
- = (internal::DataSetDescriptor<dim>::
- sub_face (cell, face, subface,
- quadrature.n_quadrature_points)).offset();
+ const typename QProjector<dim>::DataSetDescriptor offset
+ = (QProjector<dim>::DataSetDescriptor::
+ sub_face (face, subface, cell->face_orientation(face),
+ quadrature.n_quadrature_points));
const UpdateFlags flags(fe_data.update_once | fe_data.update_each);
#include <dofs/dof_accessor.h>
#include <fe/fe.h>
#include <fe/mapping.h>
-#include <fe/mapping_q1.h>
#include <fe/fe_nedelec.h>
#include <fe/fe_values.h>
if (flags & update_second_derivatives)
this->compute_2nd (mapping, cell,
- internal::DataSetDescriptor<dim>::cell().offset(),
+ QProjector<dim>::DataSetDescriptor::cell(),
mapping_data, fe_data, data);
}
// offset determines which data set
// to take (all data sets for all
// faces are stored contiguously)
- const unsigned int offset
- = (internal::DataSetDescriptor<dim>::
- face (cell, face,
- quadrature.n_quadrature_points)).offset();
+ const typename QProjector<dim>::DataSetDescriptor offset
+ = (QProjector<dim>::DataSetDescriptor::
+ face (face, cell->face_orientation(face),
+ quadrature.n_quadrature_points));
// get the flags indicating the
// fields that have to be filled
// offset determines which data set
// to take (all data sets for all
// faces are stored contiguously)
- const unsigned int offset
- = (internal::DataSetDescriptor<dim>::
- sub_face (cell, face, subface,
- quadrature.n_quadrature_points)).offset();
+ const typename QProjector<dim>::DataSetDescriptor offset
+ = (QProjector<dim>::DataSetDescriptor::
+ sub_face (face, subface, cell->face_orientation(face),
+ quadrature.n_quadrature_points));
// get the flags indicating the
// fields that have to be filled
#include <dofs/dof_accessor.h>
#include <fe/fe.h>
#include <fe/mapping.h>
-#include <fe/mapping_q1.h>
#include <fe/fe_q.h>
#include <fe/fe_values.h>
if (flags & update_second_derivatives)
this->compute_2nd (mapping, cell,
- internal::DataSetDescriptor<dim>::cell().offset(),
+ QProjector<dim>::DataSetDescriptor::cell(),
mapping_data, fe_data, data);
}
// offset determines which data set
// to take (all data sets for all
// faces are stored contiguously)
- const unsigned int offset
- = (internal::DataSetDescriptor<dim>::
- face (cell, face,
- quadrature.n_quadrature_points)).offset();
+ const typename QProjector<dim>::DataSetDescriptor offset
+ = (QProjector<dim>::DataSetDescriptor::
+ face (face, cell->face_orientation(face),
+ quadrature.n_quadrature_points));
const UpdateFlags flags(fe_data.update_once | fe_data.update_each);
// offset determines which data set
// to take (all data sets for all
// sub-faces are stored contiguously)
- const unsigned int offset
- = (internal::DataSetDescriptor<dim>::
- sub_face (cell, face, subface,
- quadrature.n_quadrature_points)).offset();
+ const typename QProjector<dim>::DataSetDescriptor offset
+ = (QProjector<dim>::DataSetDescriptor::
+ sub_face (face, subface, cell->face_orientation(face),
+ quadrature.n_quadrature_points));
const UpdateFlags flags(fe_data.update_once | fe_data.update_each);
#include <dofs/dof_accessor.h>
#include <fe/fe.h>
#include <fe/mapping.h>
-#include <fe/mapping_q1.h>
#include <fe/fe_q_hierarchical.h>
#include <fe/fe_values.h>
if (flags & update_second_derivatives)
this->compute_2nd (mapping, cell,
- internal::DataSetDescriptor<dim>::cell().offset(),
+ QProjector<dim>::DataSetDescriptor::cell(),
mapping_data, fe_data, data);
}
// offset determines which data set
// to take (all data sets for all
// faces are stored contiguously)
- const unsigned int offset
- = (internal::DataSetDescriptor<dim>::
- face (cell, face,
- quadrature.n_quadrature_points)).offset();
+ const typename QProjector<dim>::DataSetDescriptor offset
+ = (QProjector<dim>::DataSetDescriptor::
+ face (face, cell->face_orientation(face),
+ quadrature.n_quadrature_points));
const UpdateFlags flags(fe_data.update_once | fe_data.update_each);
// offset determines which data set
// to take (all data sets for all
// sub-faces are stored contiguously)
- const unsigned int offset
- = (internal::DataSetDescriptor<dim>::
- sub_face (cell, face, subface,
- quadrature.n_quadrature_points)).offset();
+ const typename QProjector<dim>::DataSetDescriptor offset
+ = (QProjector<dim>::DataSetDescriptor::
+ sub_face (face, subface, cell->face_orientation(face),
+ quadrature.n_quadrature_points));
const UpdateFlags flags(fe_data.update_once | fe_data.update_each);
#include <dofs/dof_accessor.h>
#include <fe/fe.h>
#include <fe/mapping.h>
-#include <fe/mapping_q1.h>
#include <fe/fe_raviart_thomas.h>
#include <fe/fe_values.h>
if (flags & update_second_derivatives)
this->compute_2nd (mapping, cell,
- internal::DataSetDescriptor<dim>::cell().offset(),
+ QProjector<dim>::DataSetDescriptor::cell(),
mapping_data, fe_data, data);
}
// offset determines which data set
// to take (all data sets for all
// faces are stored contiguously)
- const unsigned int offset
- = (internal::DataSetDescriptor<dim>::
- face (cell, face,
- quadrature.n_quadrature_points)).offset();
+ const typename QProjector<dim>::DataSetDescriptor offset
+ = (QProjector<dim>::DataSetDescriptor::
+ face (face, cell->face_orientation(face),
+ quadrature.n_quadrature_points));
// get the flags indicating the
// fields that have to be filled
// offset determines which data set
// to take (all data sets for all
// faces are stored contiguously)
- const unsigned int offset
- = (internal::DataSetDescriptor<dim>::
- sub_face (cell, face, subface,
- quadrature.n_quadrature_points)).offset();
+ const typename QProjector<dim>::DataSetDescriptor offset
+ = (QProjector<dim>::DataSetDescriptor::
+ sub_face (face, subface, cell->face_orientation(face),
+ quadrature.n_quadrature_points));
// get the flags indicating the
// fields that have to be filled
#include <fe/fe_raviart_thomas.h>
#include <fe/fe_system.h>
#include <fe/fe_values.h>
-#include <fe/mapping_q1.h>
#include <dofs/dof_constraints.h>
#include <dofs/dof_handler.h>
#include <dofs/dof_accessor.h>
ExcIndexRange (face_no, 0, GeometryInfo<dim>::faces_per_cell));
Assert (subface_no < GeometryInfo<dim>::subfaces_per_face,
ExcIndexRange (subface_no, 0, GeometryInfo<dim>::subfaces_per_face));
+ 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."));
this->present_cell = cell;
this->my_orientation = this->orientation_table[face_no];
#include <grid/tria_iterator.h>
#include <dofs/dof_accessor.h>
#include <fe/mapping_cartesian.h>
-#include <fe/mapping_q1.h>
#include <fe/fe_values.h>
#include <cmath>
std::vector<Point<dim> > &quadrature_points,
std::vector<Point<dim> > &normal_vectors) const
{
- UpdateFlags update_flags(data.current_update_flags());
+ const UpdateFlags update_flags(data.current_update_flags());
const unsigned int npts = quadrature_points.size ();
- typedef typename internal::DataSetDescriptor<dim> DataSetDescriptor;
- unsigned int offset = DataSetDescriptor::cell().offset();
+ const typename QProjector<dim>::DataSetDescriptor offset
+ = (face_no == invalid_face_number
+ ?
+ QProjector<dim>::DataSetDescriptor::cell()
+ :
+ (sub_no == invalid_face_number
+ ?
+ // called from FEFaceValues
+ QProjector<dim>::DataSetDescriptor::face (face_no,
+ cell->face_orientation(face_no),
+ quadrature_points.size())
+ :
+ // called from FESubfaceValues
+ QProjector<dim>::DataSetDescriptor::sub_face (face_no, sub_no,
+ cell->face_orientation(face_no),
+ quadrature_points.size())
+ ));
+ // some more sanity checks
if (face_no != invalid_face_number)
{
// Add 1 on both sides of
Assert (face_no+1 < GeometryInfo<dim>::faces_per_cell+1,
ExcIndexRange (face_no, 0, GeometryInfo<dim>::faces_per_cell));
- if (sub_no == invalid_face_number)
- // called from FEFaceValues
- offset = DataSetDescriptor::face (cell, face_no,
- quadrature_points.size()).offset();
- else
- {
- // called from
- // FESubfaceValue (do the
- // +1 trick to avoid a
- // warning about comparison
- // always being false in
- // 1d)
- Assert (sub_no+1 < GeometryInfo<dim>::subfaces_per_face+1,
- ExcIndexRange (sub_no, 0, GeometryInfo<dim>::subfaces_per_face));
- offset = DataSetDescriptor::sub_face (cell, face_no, sub_no,
- quadrature_points.size()).offset();
- }
+ Assert ((sub_no == invalid_face_number)
+ ||
+ (sub_no+1 < GeometryInfo<dim>::subfaces_per_face+1),
+ ExcIndexRange (sub_no, 0,
+ GeometryInfo<dim>::subfaces_per_face));
}
else
// invalid face number, so
J *= data.length[d];
if (data.current_update_flags() & update_JxW_values)
- {
- for (unsigned int i=0; i<JxW_values.size();++i)
- JxW_values[i] = J * q.weight(i);
- }
+ for (unsigned int i=0; i<JxW_values.size();++i)
+ JxW_values[i] = J * q.weight(i);
if (data.current_update_flags() & update_boundary_forms)
- {
- for (unsigned int i=0; i<boundary_forms.size();++i)
- boundary_forms[i] = J * normal_vectors[i];
- }
+ for (unsigned int i=0; i<boundary_forms.size();++i)
+ boundary_forms[i] = J * normal_vectors[i];
}
+
template <int dim>
void
MappingCartesian<dim>::fill_fe_subface_values (const typename DoFHandler<dim>::cell_iterator &cell,
J *= data.length[d];
if (data.current_update_flags() & update_JxW_values)
- {
- for (unsigned int i=0; i<JxW_values.size();++i)
- JxW_values[i] = J * q.weight(i) / GeometryInfo<dim>::subfaces_per_face;
- }
+ for (unsigned int i=0; i<JxW_values.size();++i)
+ JxW_values[i] = J * q.weight(i) / GeometryInfo<dim>::subfaces_per_face;
if (data.current_update_flags() & update_boundary_forms)
- {
- for (unsigned int i=0; i<boundary_forms.size();++i)
- boundary_forms[i] = J * normal_vectors[i];
- }
+ for (unsigned int i=0; i<boundary_forms.size();++i)
+ boundary_forms[i] = J * normal_vectors[i];
}
{
for (unsigned int d=0;d<dim;++d)
(*begin)[d] = (*src)[d]/data.length[d];
- begin++;
- src++;
+ ++begin;
+ ++src;
}
}
for (unsigned int d=0; d<dim; ++d)
for (unsigned int p=0; p<dim; ++p)
(*begin)[d][p] = (*src)[d][p] / data.length[p];
- begin++;
- src++;
+ ++begin;
+ ++src;
}
}
for (unsigned int d=0; d<dim; ++d)
for (unsigned int p=0; p<dim; ++p)
(*begin)[d][p] = data.length[d] * (*src)[d][p];
- begin++;
- src++;
+ ++begin;
+ ++src;
}
}
Point<dim> p_real = cell->vertex(0);
for (unsigned int d=0; d<dim; ++d)
- p_real(d) +=length[d]*p(d);
+ p_real(d) += length[d]*p(d);
return p_real;
}
#include <grid/tria_iterator.h>
#include <grid/tria_boundary.h>
#include <dofs/dof_accessor.h>
-#include <fe/mapping_q1.h>
#include <fe/fe_tools.h>
#include <numeric>
const unsigned int n_q_points=q.n_quadrature_points;
this->compute_fill_face (cell, face_no, false,
n_q_points,
- MappingQ1<dim>::DataSetDescriptor::
- face (cell, face_no, n_q_points),
+ QProjector<dim>::DataSetDescriptor::
+ face (face_no, cell->face_orientation(face_no),
+ n_q_points),
q.get_weights(),
*p_data,
quadrature_points, JxW_values,
const unsigned int n_q_points=q.n_quadrature_points;
this->compute_fill_face (cell, face_no, true,
n_q_points,
- MappingQ1<dim>::DataSetDescriptor::
- sub_face (cell, face_no, sub_no, n_q_points),
+ QProjector<dim>::DataSetDescriptor::
+ sub_face (face_no, sub_no,
+ cell->face_orientation(face_no),
+ n_q_points),
q.get_weights(),
*p_data,
quadrature_points, JxW_values,
#include <grid/tria.h>
#include <grid/tria_iterator.h>
#include <dofs/dof_accessor.h>
-#include <fe/mapping_q1.h>
#include <fe/fe_values.h>
+#include <fe/mapping_q1.h>
#include <cmath>
#include <algorithm>
-namespace internal
-{
- template <int dim>
- DataSetDescriptor<dim>
- DataSetDescriptor<dim>::cell ()
- {
- return 0;
- }
-
-
- template <int dim>
- DataSetDescriptor<dim>
- DataSetDescriptor<dim>::
- face (const typename DoFHandler<dim>::cell_iterator &cell,
- const unsigned int face_no,
- const unsigned int n_quadrature_points)
- {
- Assert (dim != 1, ExcInternalError());
- Assert (face_no < GeometryInfo<dim>::faces_per_cell,
- ExcInternalError());
- Assert (n_quadrature_points > 0, ExcInternalError());
-
- switch (dim)
- {
- case 1:
- case 2:
- return face_no * n_quadrature_points;
-
- // in 3d, we have to
- // account for faces
- // that have reverse
- // orientation. thus,
- // we have to store
- // _two_ data sets per
- // face or subface
- case 3:
- return ((face_no +
- (cell->face_orientation(face_no) == true ?
- 0 : GeometryInfo<dim>::faces_per_cell))
- * n_quadrature_points);
-
- default:
- Assert (false, ExcInternalError());
- }
- return static_cast<unsigned int>(-1);
- }
-
-
-
- template <int dim>
- DataSetDescriptor<dim>
- DataSetDescriptor<dim>::
- sub_face (const typename DoFHandler<dim>::cell_iterator &cell,
- const unsigned int face_no,
- const unsigned int subface_no,
- const unsigned int n_quadrature_points)
- {
- Assert (dim != 1, ExcInternalError());
- Assert (face_no < GeometryInfo<dim>::faces_per_cell,
- ExcInternalError());
- // the trick with +1 prevents
- // that we get a warning in 1d
- 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 1:
- case 2:
- return ((face_no * GeometryInfo<dim>::subfaces_per_face +
- subface_no)
- * n_quadrature_points);
-
- // for 3d, same as above:
- case 3:
- return (((face_no * GeometryInfo<dim>::subfaces_per_face +
- subface_no)
- + (cell->face_orientation(face_no) == true ?
- 0 :
- GeometryInfo<dim>::faces_per_cell *
- GeometryInfo<dim>::subfaces_per_face)
- )
- * n_quadrature_points);
- default:
- Assert (false, ExcInternalError());
- }
- return static_cast<unsigned int>(-1);
- }
-
-
- template <int dim>
- unsigned int
- DataSetDescriptor<dim>::offset () const
- {
- return dataset_offset;
- }
-
-
-
- template <int dim>
- DataSetDescriptor<dim>::
- DataSetDescriptor (const unsigned int dataset_offset)
- :
- dataset_offset (dataset_offset)
- {}
-
-
- template <int dim>
- DataSetDescriptor<dim>::
- DataSetDescriptor ()
- :
- dataset_offset (static_cast<unsigned int>(-1))
- {}
-}
-
-
-
template <int dim>
for (unsigned int point=0; point<n_q_points; ++point)
for (unsigned int k=0; k<data.n_shape_functions; ++k)
quadrature_points[point]
- += data.shape(point+data_set.offset(),k)
+ += data.shape(point+data_set,k)
* data.mapping_support_points[k];
// then Jacobians
for (unsigned int i=0; i<dim; ++i)
for (unsigned int j=0; j<dim; ++j)
data.contravariant[point][i][j]
- += (data.derivative(point+data_set.offset(), k)[j]
+ += (data.derivative(point+data_set, k)[j]
*
data.mapping_support_points[k][i]);
compute_fill_face (cell, face_no, false,
n_q_points,
- DataSetDescriptor::face (cell, face_no, n_q_points),
+ DataSetDescriptor::face (face_no,
+ cell->face_orientation(face_no),
+ n_q_points),
q.get_weights(),
data,
quadrature_points,
compute_fill_face (cell, face_no, true,
n_q_points,
- DataSetDescriptor::sub_face (cell, face_no, sub_no, n_q_points),
+ DataSetDescriptor::sub_face (face_no, sub_no,
+ cell->face_orientation(face_no),
+ n_q_points),
q.get_weights(),
data,
quadrature_points,
//----------------------------------------------------------------------//
template class MappingQ1<deal_II_dimension>;
-namespace internal
-{
- template class DataSetDescriptor<deal_II_dimension>;
-}
const active_cell_iterator neighbor = cell->neighbor(face_no);
- // find which number the current
- // face has relative to the neighboring
- // cell
- const unsigned int neighbor_neighbor = cell->neighbor_of_neighbor (face_no);
- Assert (neighbor_neighbor<GeometryInfo<dim>::faces_per_cell, ExcInternalError());
+ // find which number the
+ // current face has relative to
+ // the neighboring cell
+ const unsigned int neighbor_neighbor
+ = cell->neighbor_of_neighbor (face_no);
+ Assert (neighbor_neighbor<GeometryInfo<dim>::faces_per_cell,
+ ExcInternalError());
// get restriction of finite element
// function of @p{neighbor} to the
// index the number of the
// quadrature point
- // store which number @p{cell} has in the
- // list of neighbors of @p{neighbor}
- const unsigned int neighbor_neighbor = cell->neighbor_of_neighbor (face_no);
- Assert (neighbor_neighbor<GeometryInfo<dim>::faces_per_cell, ExcInternalError());
+ // store which number @p{cell} has
+ // in the list of neighbors of
+ // @p{neighbor}
+ const unsigned int neighbor_neighbor
+ = cell->neighbor_of_neighbor (face_no);
+ Assert (neighbor_neighbor<GeometryInfo<dim>::faces_per_cell,
+ ExcInternalError());
// loop over all subfaces
- for (unsigned int subface_no=0; subface_no<GeometryInfo<dim>::subfaces_per_face;
+ for (unsigned int subface_no=0;
+ subface_no<GeometryInfo<dim>::subfaces_per_face;
++subface_no)
{
// get an iterator pointing to the
// cell behind the present subface
- 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,
- (face_orientations_match ?
- subface_no :
- subface_translation[subface_no])));
const active_cell_iterator neighbor_child
- = neighbor->child(neighbor_child_index);
- Assert (!neighbor->child(neighbor_child_index)->has_children(),
+ = cell->neighbor_child_on_subface (face_no, subface_no);
+ Assert (!neighbor_child->has_children(),
ExcInternalError());
// restrict the finite element
<h3 style="color:red">Incompatibilities</h3>
<ol>
+ <li> <p>
+ Changed: The <code class="class">QProjector</class> has functions that
+ project a lower-dimensional quadrature formula onto all faces or
+ subfaces of a cell. In 3d, it now does this but also adds projections of
+ these quadrature formula onto the faces from the other side. We need
+ this due to the fact that we now support faces in 3d that have a normal
+ vector opposite to the standard direction.
+ <br>
+ (WB 2003/10/17)
+ </p>
+
<li> <p> Moved and changed: The header file
- <tt>include/numerics/dof_renumbering.h</tt> has been moved to the
- directory <tt>include/dofs</tt>, where it logically
- belongs. Furthermore, the sorting parameter of the function <code
- class="class">DoFRenumbering</code>::<code
- class="member">component_wise</code> has changed its meaning. See
- the reference documentation for details.
- <br>
- (GK 2003/07/03)
+ <tt>include/numerics/dof_renumbering.h</tt> has been moved to the
+ directory <tt>include/dofs</tt>, where it logically
+ belongs. Furthermore, the sorting parameter of the function <code
+ class="class">DoFRenumbering</code>::<code
+ class="member">component_wise</code> has changed its meaning. See
+ the reference documentation for details.
+ <br>
+ (GK 2003/07/03)
</ol>
<h3>General</h3>
<ol>
+ <li> <p>
+ Augmented: The <code
+ class="member">GeometryInfo::child_cell_on_face</code>
+ results in a result that might not be what you expect in 3d in some
+ cases. A warning has been added to the documentation, as well as a
+ reference to the function to call instead.
+ <br>
+ (WB 2003/10/17)
+ </p>
+
<li> <p> Fixed: The step-14 program had a bug in the rare case that
there are more CPUs in a machine than there are cells. This is
now fixed.
<h3>deal.II</h3>
<ol>
+ <li> <p>
+ New: There is now a function <code
+ class="member">CellAccessor::neighbor_child_on_subface</code>
+ that should be called instead of using <code
+ class="member">GeometryInfo::child_cell_on_face</code> in most cases.
+ <br>
+ (WB 2003/10/17)
+ </p>
+
<li> <p>
New: <code class="class">GridGenerator</code> has a new
function <code class="member">subdivided_hyper_rectangle</code>
// `behind' the
// current
// subface.
- typename DoFHandler<dim>::active_cell_iterator neighbor_child=
- neighbor->child(GeometryInfo<dim>::
- child_cell_on_face(neighbor2,subface_no));
+ typename DoFHandler<dim>::active_cell_iterator
+ neighbor_child
+ = cell->neighbor_child_on_subface (face_no, subface_no);
// As these are
// quite
const unsigned int neighbor_face_no=faceno_subfaceno.first,
neighbor_subface_no=faceno_subfaceno.second;
- Assert (neighbor->neighbor(neighbor_face_no)
- ->child(GeometryInfo<dim>::child_cell_on_face(
- face_no,neighbor_subface_no)) == cell, ExcInternalError());
+ Assert (neighbor->neighbor_child_on_subface (neighbor_face_no,
+ neighbor_subface_no)
+ == cell,
+ ExcInternalError());
// Reinit the
// appropriate
subface_no<GeometryInfo<dim>::subfaces_per_face;
++subface_no)
{
- typename DoFHandler<dim>::cell_iterator neighbor_child=
- neighbor->child(GeometryInfo<dim>::child_cell_on_face(
- neighbor2,subface_no));
+ typename DoFHandler<dim>::cell_iterator neighbor_child
+ = cell->neighbor_child_on_subface (face_no, subface_no);
Assert (neighbor_child->face(neighbor2) == face->child(subface_no),
ExcInternalError());
Assert (!neighbor_child->has_children(), ExcInternalError());
// assertion will be removed
// anyway.
const active_cell_iterator neighbor_child
- = neighbor->child(GeometryInfo<dim>::
- child_cell_on_face(neighbor_neighbor,
- subface_no));
+ = cell->neighbor_child_on_subface (face_no, subface_no);
Assert (neighbor_child->face(neighbor_neighbor) ==
cell->face(face_no)->child(subface_no),
ExcInternalError());