get_prolongation_matrix (const unsigned int child,
const RefinementCase<dim> &refinement_case=RefinementCase<dim>::isotropic_refinement) const;
+ /**
+ * Given an index in the natural ordering of indices on a face, return the
+ * index of the same degree of freedom on the cell.
+ *
+ * To explain the concept, consider the case where we would like to know
+ * whether a degree of freedom on a face, for example as part of an FESystem
+ * element, is primitive. Unfortunately, the
+ * is_primitive() function in the FiniteElement class takes a cell index, so
+ * we would need to find the cell index of the shape function that
+ * corresponds to the present face index. This function does that.
+ *
+ * Code implementing this would then look like this:
+ * @code
+ * for (i=0; i<dofs_per_face; ++i)
+ * if (fe.is_primitive(fe.face_to_equivalent_cell_index(i, some_face_no)))
+ * ... do whatever
+ * @endcode
+ * The function takes additional arguments that account for the fact that
+ * actual faces can be in their standard ordering with respect to the cell
+ * under consideration, or can be flipped, oriented, etc.
+ *
+ * @param face_dof_index The index of the degree of freedom on a face.
+ * This index must be between zero and dofs_per_face.
+ * @param face The number of the face this degree of freedom lives on.
+ * This number must be between zero and GeometryInfo::faces_per_cell.
+ * @param face_orientation One part of the description of the orientation
+ * of the face. See @ref GlossFaceOrientation .
+ * @param face_flip One part of the description of the orientation
+ * of the face. See @ref GlossFaceOrientation .
+ * @param face_rotation One part of the description of the orientation
+ * of the face. See @ref GlossFaceOrientation .
+ * @return The index of this degree of freedom within the set
+ * of degrees of freedom on the entire cell. The returned value
+ * will be between zero and dofs_per_cell.
+ */
+ virtual
+ unsigned int face_to_cell_index (const unsigned int face_dof_index,
+ const unsigned int face,
+ const bool face_orientation = true,
+ const bool face_flip = false,
+ const bool face_rotation = false) const;
+
/**
* @name Functions to support hp
* @{
"an overloaded version but apparently hasn't done so. See "
"the documentation of this function for more information."));
- // DoF on a vertex
+ // we need to distinguish between DoFs on vertices, lines and in 3d quads.
+ // do so in a sequence of if-else statements
if (face_index < this->first_face_line_index)
+ // DoF is on a vertex
{
- // Vertex number on the face
- const unsigned int face_vertex = face_index / this->dofs_per_vertex;
- return face_index % this->dofs_per_vertex
- + GeometryInfo<dim>::face_to_cell_vertices(face, face_vertex,
- face_orientation,
- face_flip,
- face_rotation)
- * this->dofs_per_vertex;
+ // get the number of the vertex on the face that corresponds to this DoF,
+ // along with the number of the DoF on this vertex
+ const unsigned int face_vertex = face_index / this->dofs_per_vertex;
+ const unsigned int dof_index_on_vertex = face_index % this->dofs_per_vertex;
+
+ // then get the number of this vertex on the cell and translate
+ // this to a DoF number on the cell
+ return (GeometryInfo<dim>::face_to_cell_vertices(face, face_vertex,
+ face_orientation,
+ face_flip,
+ face_rotation)
+ * this->dofs_per_vertex
+ +
+ dof_index_on_vertex);
}
- // Else, DoF on a line?
- if (face_index < this->first_face_quad_index)
+ else if (face_index < this->first_face_quad_index)
+ // DoF is on a face
{
- // Ignore vertex dofs
+ // do the same kind of translation as before. we need to only consider
+ // DoFs on the lines, i.e., ignoring those on the vertices
const unsigned int index = face_index - this->first_face_line_index;
- // Line number on the face
- const unsigned int face_line = index / this->dofs_per_line;
- return this->first_line_index + index % this->dofs_per_line
- + GeometryInfo<dim>::face_to_cell_lines(face, face_line,
- face_orientation,
- face_flip,
- face_rotation)
- * this->dofs_per_line;
+
+ const unsigned int face_line = index / this->dofs_per_line;
+ const unsigned int dof_index_on_line = index % this->dofs_per_line;
+
+ return (this->first_line_index
+ + GeometryInfo<dim>::face_to_cell_lines(face, face_line,
+ face_orientation,
+ face_flip,
+ face_rotation)
+ * this->dofs_per_line
+ +
+ dof_index_on_line);
}
- // Else, DoF is on a quad
+ else
+ // DoF is on a quad
+ {
+ Assert (dim >= 3, ExcInternalError());
- // Ignore vertex and line dofs
- const unsigned int index = face_index - this->first_face_quad_index;
- return this->first_quad_index + index
- + face * this->dofs_per_quad;
+ // ignore vertex and line dofs
+ const unsigned int index = face_index - this->first_face_quad_index;
+
+ return (this->first_quad_index
+ + face * this->dofs_per_quad
+ + index);
+ }
}
+template <class POLY, int dim, int spacedim>
+unsigned int
+FE_Q_Base<POLY,dim,spacedim>::
+face_to_cell_index (const unsigned int face_index,
+ const unsigned int face,
+ const bool face_orientation,
+ const bool face_flip,
+ const bool face_rotation) const
+{
+ Assert (face_index < this->dofs_per_face,
+ ExcIndexRange(face_index, 0, this->dofs_per_face));
+ Assert (face < GeometryInfo<dim>::faces_per_cell,
+ ExcIndexRange(face, 0, GeometryInfo<dim>::faces_per_cell));
+
+ // we need to distinguish between DoFs on vertices, lines and in 3d quads.
+ // do so in a sequence of if-else statements
+ if (face_index < this->first_face_line_index)
+ // DoF is on a vertex
+ {
+ // get the number of the vertex on the face that corresponds to this DoF,
+ // along with the number of the DoF on this vertex
+ const unsigned int face_vertex = face_index / this->dofs_per_vertex;
+ const unsigned int dof_index_on_vertex = face_index % this->dofs_per_vertex;
+
+ // then get the number of this vertex on the cell and translate
+ // this to a DoF number on the cell
+ return (GeometryInfo<dim>::face_to_cell_vertices(face, face_vertex,
+ face_orientation,
+ face_flip,
+ face_rotation)
+ * this->dofs_per_vertex
+ +
+ dof_index_on_vertex);
+ }
+ else if (face_index < this->first_face_quad_index)
+ // DoF is on a face
+ {
+ // do the same kind of translation as before. we need to only consider
+ // DoFs on the lines, i.e., ignoring those on the vertices
+ const unsigned int index = face_index - this->first_face_line_index;
+
+ const unsigned int face_line = index / this->dofs_per_line;
+ const unsigned int dof_index_on_line = index % this->dofs_per_line;
+
+ // we now also need to adjust the line index for the case of
+ // face orientation, face flips, etc
+ unsigned int adjusted_dof_index_on_line;
+ switch (dim)
+ {
+ case 1:
+ Assert (false, ExcInternalError());
+
+ case 2:
+ // in 2d, only face_flip has a meaning. if it is set, consider
+ // dofs in reverse order
+ if (face_flip == false)
+ adjusted_dof_index_on_line = dof_index_on_line;
+ else
+ adjusted_dof_index_on_line = this->dofs_per_line - dof_index_on_line - 1;
+ break;
+
+ case 3:
+ // in 3d, things are difficult. someone will have to think
+ // about how this code here should look like, by drawing a bunch
+ // of pictures of how all the faces can look like with the various
+ // flips and rotations.
+ //
+ // that said, the Q2 case is easy enough to implement
+ Assert (this->dofs_per_line <= 1, ExcNotImplemented());
+ adjusted_dof_index_on_line = dof_index_on_line;
+ break;
+ }
+
+ return (this->first_line_index
+ + GeometryInfo<dim>::face_to_cell_lines(face, face_line,
+ face_orientation,
+ face_flip,
+ face_rotation)
+ * this->dofs_per_line
+ +
+ adjusted_dof_index_on_line);
+ }
+ else
+ // DoF is on a quad
+ {
+ Assert (dim >= 3, ExcInternalError());
+
+ // ignore vertex and line dofs
+ const unsigned int index = face_index - this->first_face_quad_index;
+
+ // the same is true here as above for the 3d case -- someone will
+ // just have to draw a bunch of pictures. in the meantime,
+ // we can implement the Q2 case in which it is simple
+ Assert (this->dofs_per_quad <= 1, ExcNotImplemented());
+ return (this->first_quad_index
+ + face * this->dofs_per_quad
+ + index);
+ }
+}
+
+
+
template <class POLY, int dim, int spacedim>
std::vector<unsigned int>
--- /dev/null
+// ---------------------------------------------------------------------
+// $Id$
+//
+// Copyright (C) 1998 - 2013 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+
+// it turns out that FE_Q::face_to_cell_index() had a bug for elements beyond
+// Q2 when using the face flip flag. this test is for the 2d case for the Q3
+// case
+
+#include "../tests.h"
+#include <iostream>
+#include <fstream>
+
+#include <deal.II/base/logstream.h>
+#include <deal.II/fe/fe_q.h>
+
+template<int dim>
+void test()
+{
+ FE_Q<dim> fe(3);
+ const unsigned int dofs_per_face = fe.dofs_per_face;
+
+ for (unsigned int face=0; face<4; ++face)
+ {
+ deallog << "Face=" << face << std::endl;
+
+ for (int flip=0; flip<2; ++flip)
+ {
+ deallog << " flip=" << (flip == 0 ? "false" : "true")
+ << std::endl
+ << " ";
+ for (unsigned int i = 0; i < dofs_per_face; ++i)
+ deallog << fe.face_to_cell_index(i, face,
+ true,
+ (flip == 0 ? false : true),
+ false) << " - ";
+ deallog << std::endl;
+ }
+ }
+}
+
+int main()
+{
+ std::ofstream logfile("face_to_cell_q3_2d/output");
+ deallog.attach(logfile);
+ deallog.depth_console(0);
+ deallog.threshold_double(1.e-10);
+
+ test<2>();
+}
+
+
+
--- /dev/null
+
+DEAL::Face=0
+DEAL:: flip=false
+DEAL:: 0 - 2 - 4 - 5 -
+DEAL:: flip=true
+DEAL:: 2 - 0 - 5 - 4 -
+DEAL::Face=1
+DEAL:: flip=false
+DEAL:: 1 - 3 - 6 - 7 -
+DEAL:: flip=true
+DEAL:: 3 - 1 - 7 - 6 -
+DEAL::Face=2
+DEAL:: flip=false
+DEAL:: 0 - 1 - 8 - 9 -
+DEAL:: flip=true
+DEAL:: 1 - 0 - 9 - 8 -
+DEAL::Face=3
+DEAL:: flip=false
+DEAL:: 2 - 3 - 10 - 11 -
+DEAL:: flip=true
+DEAL:: 3 - 2 - 11 - 10 -
--- /dev/null
+// ---------------------------------------------------------------------
+// $Id$
+//
+// Copyright (C) 1998 - 2013 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+
+// it turns out that FE_Q::face_to_cell_index() had a bug for elements beyond
+// Q2 when using the face flip flag. this test is for the 2d case for the Q4
+// case
+
+#include "../tests.h"
+#include <iostream>
+#include <fstream>
+
+#include <deal.II/base/logstream.h>
+#include <deal.II/fe/fe_q.h>
+
+template<int dim>
+void test()
+{
+ FE_Q<dim> fe(4);
+ const unsigned int dofs_per_face = fe.dofs_per_face;
+
+ for (unsigned int face=0; face<4; ++face)
+ {
+ deallog << "Face=" << face << std::endl;
+
+ for (int flip=0; flip<2; ++flip)
+ {
+ deallog << " flip=" << (flip == 0 ? "false" : "true")
+ << std::endl
+ << " ";
+ for (unsigned int i = 0; i < dofs_per_face; ++i)
+ deallog << fe.face_to_cell_index(i, face,
+ true,
+ (flip == 0 ? false : true),
+ false) << " - ";
+ deallog << std::endl;
+ }
+ }
+}
+
+int main()
+{
+ std::ofstream logfile("face_to_cell_q4_2d/output");
+ deallog.attach(logfile);
+ deallog.depth_console(0);
+ deallog.threshold_double(1.e-10);
+
+ test<2>();
+}
+
+
+
--- /dev/null
+
+DEAL::Face=0
+DEAL:: flip=false
+DEAL:: 0 - 2 - 4 - 5 - 6 -
+DEAL:: flip=true
+DEAL:: 2 - 0 - 6 - 5 - 4 -
+DEAL::Face=1
+DEAL:: flip=false
+DEAL:: 1 - 3 - 7 - 8 - 9 -
+DEAL:: flip=true
+DEAL:: 3 - 1 - 9 - 8 - 7 -
+DEAL::Face=2
+DEAL:: flip=false
+DEAL:: 0 - 1 - 10 - 11 - 12 -
+DEAL:: flip=true
+DEAL:: 1 - 0 - 12 - 11 - 10 -
+DEAL::Face=3
+DEAL:: flip=false
+DEAL:: 2 - 3 - 13 - 14 - 15 -
+DEAL:: flip=true
+DEAL:: 3 - 2 - 15 - 14 - 13 -