]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Implement face_to_cell_index() also for FESystem.
authorWolfgang Bangerth <bangerth@math.tamu.edu>
Sat, 10 Aug 2013 11:44:30 +0000 (11:44 +0000)
committerWolfgang Bangerth <bangerth@math.tamu.edu>
Sat, 10 Aug 2013 11:44:30 +0000 (11:44 +0000)
git-svn-id: https://svn.dealii.org/trunk@30277 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/doc/news/changes.h
deal.II/include/deal.II/fe/fe_system.h
deal.II/source/fe/fe_system.cc
tests/fe/face_to_cell_q2xq2_2d.cc [new file with mode: 0644]
tests/fe/face_to_cell_q2xq2_2d/cmp/generic [new file with mode: 0644]
tests/fe/face_to_cell_q3xq4_2d.cc [new file with mode: 0644]
tests/fe/face_to_cell_q3xq4_2d/cmp/generic [new file with mode: 0644]

index ad91e1e9855951466bd97d309c98d3a465761bc6..b02cded52bb3e0c1fd6f178802035a34f2920aa5 100644 (file)
@@ -1,5 +1,5 @@
 /**
- * @page changes_after_8_0 Changes after Version 8.0
+// * @page changes_after_8_0 Changes after Version 8.0
 
 <p>
 This is the list of changes made after the release of
@@ -54,6 +54,18 @@ back already. Please use FiniteElement::face_to_cell_index() instead.
 <h3>Specific improvements</h3>
 
 <ol>
+  <li>
+  Fixed: The FiniteElement::face_to_cell_index() function had a bug
+  that made it work incorrectly for elements that have more than one
+  degree of freedom per line (in 2d) or per quad (in 3d). This is now
+  fixed for the most common cases, namely the FE_Q elements as well
+  as elements composed of FESystem elements. For all other cases, an
+  exception is generated reporting that this case is not implemented.
+  If you run into this, let us know.
+  <br>
+  (Wolfgang Bangerth, 2013/08/10)
+  </li>
+
   <li>
   New: DataOutBase::VtkFlags now has a flag
   DataOutBase::VtkFlags::print_date_and_time that can be used to suppress output
index 96a0d2dfbd80c2ad86d849a83b1ed751e273307d..083a167da22d12d603c1929db7e3730517da3752 100644 (file)
@@ -413,6 +413,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;
+
   /**
    * Implementation of the respective function in the base class.
    */
index 2f45064642d9fd443280e4f991af1864fedddb1f..3df96eddb335b7abf70cdcce6590de806e332ca5 100644 (file)
@@ -745,6 +745,45 @@ FESystem<dim,spacedim>
 }
 
 
+template <int dim, int spacedim>
+unsigned int
+FESystem<dim,spacedim>::
+face_to_cell_index (const unsigned int face_dof_index,
+                    const unsigned int face,
+                    const bool face_orientation,
+                    const bool face_flip,
+                    const bool face_rotation) const
+{
+  // we need to ask the base elements how they want to translate
+  // the DoFs within their own numbering. thus, translate to
+  // the base element numbering and then back
+  const std::pair<std::pair<unsigned int, unsigned int>, unsigned int>
+  face_base_index = this->face_system_to_base_index(face_dof_index);
+
+  const unsigned int
+  base_face_to_cell_index
+  = this->base_element(face_base_index.first.first).face_to_cell_index (face_base_index.second,
+                                                                        face,
+                                                                        face_orientation,
+                                                                        face_flip,
+                                                                        face_rotation);
+
+  // it would be nice if we had a base_to_system_index function, but
+  // all that exists is a component_to_system_index function. we can't do
+  // this here because it won't work for non-primitive elements. consequently,
+  // simply do a loop over all dofs till we find whether it corresponds
+  // to the one we're interested in -- crude, maybe, but works for now
+  const std::pair<std::pair<unsigned int, unsigned int>, unsigned int>
+  target = std::make_pair (face_base_index.first, base_face_to_cell_index);
+  for (unsigned int i=0; i<this->dofs_per_cell; ++i)
+    if (this->system_to_base_index(i) == target)
+      return i;
+
+  Assert (false, ExcInternalError());
+  return numbers::invalid_unsigned_int;
+}
+
+
 
 //---------------------------------------------------------------------------
 // Data field initialization
diff --git a/tests/fe/face_to_cell_q2xq2_2d.cc b/tests/fe/face_to_cell_q2xq2_2d.cc
new file mode 100644 (file)
index 0000000..a515f9a
--- /dev/null
@@ -0,0 +1,67 @@
+// ---------------------------------------------------------------------
+// $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
+// Q2XQ2 when using the face flip flag. this test is for the 2d case for the Q2XQ2
+// case
+
+#include "../tests.h"
+#include <iostream>
+#include <fstream>
+
+#include <deal.II/base/logstream.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_system.h>
+
+template<int dim>
+void test()
+{
+  FESystem<dim> fe(FE_Q<dim>(2), 1,
+                  FE_Q<dim>(2), 1);
+  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_q2xq2_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_q2xq2_2d/cmp/generic b/tests/fe/face_to_cell_q2xq2_2d/cmp/generic
new file mode 100644 (file)
index 0000000..15b96ad
--- /dev/null
@@ -0,0 +1,21 @@
+
+DEAL::Face=0
+DEAL::  flip=false
+DEAL::    0 - 1 - 4 - 5 - 8 - 9 - 
+DEAL::  flip=true
+DEAL::    4 - 5 - 0 - 1 - 8 - 9 - 
+DEAL::Face=1
+DEAL::  flip=false
+DEAL::    2 - 3 - 6 - 7 - 10 - 11 - 
+DEAL::  flip=true
+DEAL::    6 - 7 - 2 - 3 - 10 - 11 - 
+DEAL::Face=2
+DEAL::  flip=false
+DEAL::    0 - 1 - 2 - 3 - 12 - 13 - 
+DEAL::  flip=true
+DEAL::    2 - 3 - 0 - 1 - 12 - 13 - 
+DEAL::Face=3
+DEAL::  flip=false
+DEAL::    4 - 5 - 6 - 7 - 14 - 15 - 
+DEAL::  flip=true
+DEAL::    6 - 7 - 4 - 5 - 14 - 15 - 
diff --git a/tests/fe/face_to_cell_q3xq4_2d.cc b/tests/fe/face_to_cell_q3xq4_2d.cc
new file mode 100644 (file)
index 0000000..52c059e
--- /dev/null
@@ -0,0 +1,67 @@
+// ---------------------------------------------------------------------
+// $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 Q3XQ4
+// case
+
+#include "../tests.h"
+#include <iostream>
+#include <fstream>
+
+#include <deal.II/base/logstream.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_system.h>
+
+template<int dim>
+void test()
+{
+  FESystem<dim> fe(FE_Q<dim>(3), 1,
+                  FE_Q<dim>(4), 1);
+  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_q3xq4_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_q3xq4_2d/cmp/generic b/tests/fe/face_to_cell_q3xq4_2d/cmp/generic
new file mode 100644 (file)
index 0000000..3f20c26
--- /dev/null
@@ -0,0 +1,21 @@
+
+DEAL::Face=0
+DEAL::  flip=false
+DEAL::    0 - 1 - 4 - 5 - 8 - 9 - 10 - 11 - 12 - 
+DEAL::  flip=true
+DEAL::    4 - 5 - 0 - 1 - 9 - 8 - 12 - 11 - 10 - 
+DEAL::Face=1
+DEAL::  flip=false
+DEAL::    2 - 3 - 6 - 7 - 13 - 14 - 15 - 16 - 17 - 
+DEAL::  flip=true
+DEAL::    6 - 7 - 2 - 3 - 14 - 13 - 17 - 16 - 15 - 
+DEAL::Face=2
+DEAL::  flip=false
+DEAL::    0 - 1 - 2 - 3 - 18 - 19 - 20 - 21 - 22 - 
+DEAL::  flip=true
+DEAL::    2 - 3 - 0 - 1 - 19 - 18 - 22 - 21 - 20 - 
+DEAL::Face=3
+DEAL::  flip=false
+DEAL::    4 - 5 - 6 - 7 - 23 - 24 - 25 - 26 - 27 - 
+DEAL::  flip=true
+DEAL::    6 - 7 - 4 - 5 - 24 - 23 - 27 - 26 - 25 - 

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.