From: bangerth Date: Sat, 10 Aug 2013 10:52:57 +0000 (+0000) Subject: Implement the face_to_cell() function for the FE_Q element family, at least for 2d... X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=d1bbf94badfc57ba16def19a2d2537a447c1e9ab;p=dealii-svn.git Implement the face_to_cell() function for the FE_Q element family, at least for 2d. In 3d, someone will need to draw a few pictures and figure things out, but at least for now one gets an ExcNotImplemented(). git-svn-id: https://svn.dealii.org/trunk@30276 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/include/deal.II/fe/fe_q_base.h b/deal.II/include/deal.II/fe/fe_q_base.h index 0bbeb3bf73..15b6701171 100644 --- a/deal.II/include/deal.II/fe/fe_q_base.h +++ b/deal.II/include/deal.II/fe/fe_q_base.h @@ -153,6 +153,48 @@ public: get_prolongation_matrix (const unsigned int child, const RefinementCase &refinement_case=RefinementCase::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; ifirst_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::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::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::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::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); + } } diff --git a/deal.II/source/fe/fe_q_base.cc b/deal.II/source/fe/fe_q_base.cc index 71b776635c..36491b4237 100644 --- a/deal.II/source/fe/fe_q_base.cc +++ b/deal.II/source/fe/fe_q_base.cc @@ -986,6 +986,108 @@ FE_Q_Base::initialize_quad_dof_index_permutation () +template +unsigned int +FE_Q_Base:: +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::faces_per_cell, + ExcIndexRange(face, 0, GeometryInfo::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::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::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 std::vector diff --git a/tests/fe/face_to_cell_q3_2d.cc b/tests/fe/face_to_cell_q3_2d.cc new file mode 100644 index 0000000000..cbbb5ffe75 --- /dev/null +++ b/tests/fe/face_to_cell_q3_2d.cc @@ -0,0 +1,65 @@ +// --------------------------------------------------------------------- +// $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 +#include + +#include +#include + +template +void test() +{ + FE_Q 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>(); +} + + + diff --git a/tests/fe/face_to_cell_q3_2d/cmp/generic b/tests/fe/face_to_cell_q3_2d/cmp/generic new file mode 100644 index 0000000000..6cac1f4c4d --- /dev/null +++ b/tests/fe/face_to_cell_q3_2d/cmp/generic @@ -0,0 +1,21 @@ + +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 - diff --git a/tests/fe/face_to_cell_q4_2d.cc b/tests/fe/face_to_cell_q4_2d.cc new file mode 100644 index 0000000000..f1e04c0b2f --- /dev/null +++ b/tests/fe/face_to_cell_q4_2d.cc @@ -0,0 +1,65 @@ +// --------------------------------------------------------------------- +// $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 +#include + +#include +#include + +template +void test() +{ + FE_Q 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>(); +} + + + diff --git a/tests/fe/face_to_cell_q4_2d/cmp/generic b/tests/fe/face_to_cell_q4_2d/cmp/generic new file mode 100644 index 0000000000..81ffa7c977 --- /dev/null +++ b/tests/fe/face_to_cell_q4_2d/cmp/generic @@ -0,0 +1,21 @@ + +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 -