]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Implement the face_to_cell() function for the FE_Q element family, at least for 2d...
authorWolfgang Bangerth <bangerth@math.tamu.edu>
Sat, 10 Aug 2013 10:52:57 +0000 (10:52 +0000)
committerWolfgang Bangerth <bangerth@math.tamu.edu>
Sat, 10 Aug 2013 10:52:57 +0000 (10:52 +0000)
git-svn-id: https://svn.dealii.org/trunk@30276 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/include/deal.II/fe/fe_q_base.h
deal.II/source/fe/fe.cc
deal.II/source/fe/fe_q_base.cc
tests/fe/face_to_cell_q3_2d.cc [new file with mode: 0644]
tests/fe/face_to_cell_q3_2d/cmp/generic [new file with mode: 0644]
tests/fe/face_to_cell_q4_2d.cc [new file with mode: 0644]
tests/fe/face_to_cell_q4_2d/cmp/generic [new file with mode: 0644]

index 0bbeb3bf738b5f19146b0bf16bc390429c1a1c02..15b6701171bd402597522306c8dc59159b7aa36a 100644 (file)
@@ -153,6 +153,48 @@ public:
   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
    * @{
index f5332eb778801ff6e158a983eb468695e1a5edf0..8e48a01f87df2dcce1cb29a325962c9fb9f0e3e4 100644 (file)
@@ -569,38 +569,57 @@ face_to_cell_index (const unsigned int face_index,
                         "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);
+    }
 }
 
 
index 71b776635c5ee816698e5a1293f4677b1a46f363..36491b4237b4592e8de6e3a14a2abd957fb63b1d 100644 (file)
@@ -986,6 +986,108 @@ FE_Q_Base<POLY,dim,spacedim>::initialize_quad_dof_index_permutation ()
 
 
 
+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>
diff --git a/tests/fe/face_to_cell_q3_2d.cc b/tests/fe/face_to_cell_q3_2d.cc
new file mode 100644 (file)
index 0000000..cbbb5ff
--- /dev/null
@@ -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 <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>();
+}
+
+
+
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 (file)
index 0000000..6cac1f4
--- /dev/null
@@ -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 (file)
index 0000000..f1e04c0
--- /dev/null
@@ -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 <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>();
+}
+
+
+
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 (file)
index 0000000..81ffa7c
--- /dev/null
@@ -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 - 

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.