]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Allow Q(p) and DGQ(r) elements to touch each other in hp::DoFHandlers. 175/head
authorWolfgang Bangerth <bangerth@math.tamu.edu>
Fri, 26 Sep 2014 15:21:10 +0000 (10:21 -0500)
committerWolfgang Bangerth <bangerth@math.tamu.edu>
Sat, 27 Sep 2014 11:38:53 +0000 (06:38 -0500)
The thing is that the FE_DGQ is discontinuous anyway, so there can't be any constraints to begin with.

doc/news/changes.h
source/fe/fe_dgq.cc
source/fe/fe_q_base.cc
tests/hp/q_vs_dgq.cc [new file with mode: 0644]
tests/hp/q_vs_dgq.output [new file with mode: 0644]

index b3829b40abdb2cfd61243bb4b8e627786531ba99..d7d5a6b361f91ec5bd184b59b85b266a2de253d2 100644 (file)
@@ -295,6 +295,15 @@ inconvenience this causes.
 <h3>Specific improvements</h3>
 
 <ol>
+  <li> Fixed: Trying to have FE_Q(p) and FE_DGQ(r) elements next to each
+  other in an hp::DoFHandler object led to assertions saying that these two
+  elements don't know how to compute interface constraints where such
+  elements touch. This has now been fixed: FE_DGQ is a discontinuous element,
+  so there cannot be any interface constraints at all.
+  <br>
+  (Wolfgang Bangerth, 2014/09/26)
+  </li>
+
   <li> Fixed: The function TrilinosWrappers::VectorBase::sadd(double factor,
   VectorBase &v) erroneously added factor*v instead of scaling the calling
   vector by factor. This is now fixed.
index a61ecfb8ae3000e6fe02796bd67565b47c5f7f92..b00486470a2d13e59892eda961b5d874c434b52b 100644 (file)
@@ -617,16 +617,12 @@ std::vector<std::pair<unsigned int, unsigned int> >
 FE_DGQ<dim, spacedim>::
 hp_vertex_dof_identities (const FiniteElement<dim, spacedim> &fe_other) const
 {
-  // there are no such constraints for DGQ
-  // elements at all
-  if (dynamic_cast<const FE_DGQ<dim, spacedim>*>(&fe_other) != 0)
-    return
-      std::vector<std::pair<unsigned int, unsigned int> > ();
-  else
-    {
-      Assert (false, ExcNotImplemented());
-      return std::vector<std::pair<unsigned int, unsigned int> > ();
-    }
+  // this element is discontinuous, so by definition there can
+  // be no identities between its dofs and those of any neighbor
+  // (of whichever type the neighbor may be -- after all, we have
+  // no face dofs on this side to begin with)
+  return
+    std::vector<std::pair<unsigned int, unsigned int> > ();
 }
 
 
@@ -636,16 +632,12 @@ std::vector<std::pair<unsigned int, unsigned int> >
 FE_DGQ<dim, spacedim>::
 hp_line_dof_identities (const FiniteElement<dim, spacedim> &fe_other) const
 {
-  // there are no such constraints for DGQ
-  // elements at all
-  if (dynamic_cast<const FE_DGQ<dim, spacedim>*>(&fe_other) != 0)
-    return
-      std::vector<std::pair<unsigned int, unsigned int> > ();
-  else
-    {
-      Assert (false, ExcNotImplemented());
-      return std::vector<std::pair<unsigned int, unsigned int> > ();
-    }
+  // this element is discontinuous, so by definition there can
+  // be no identities between its dofs and those of any neighbor
+  // (of whichever type the neighbor may be -- after all, we have
+  // no face dofs on this side to begin with)
+  return
+    std::vector<std::pair<unsigned int, unsigned int> > ();
 }
 
 
@@ -655,16 +647,12 @@ std::vector<std::pair<unsigned int, unsigned int> >
 FE_DGQ<dim, spacedim>::
 hp_quad_dof_identities (const FiniteElement<dim, spacedim>        &fe_other) const
 {
-  // there are no such constraints for DGQ
-  // elements at all
-  if (dynamic_cast<const FE_DGQ<dim, spacedim>*>(&fe_other) != 0)
-    return
-      std::vector<std::pair<unsigned int, unsigned int> > ();
-  else
-    {
-      Assert (false, ExcNotImplemented());
-      return std::vector<std::pair<unsigned int, unsigned int> > ();
-    }
+  // this element is discontinuous, so by definition there can
+  // be no identities between its dofs and those of any neighbor
+  // (of whichever type the neighbor may be -- after all, we have
+  // no face dofs on this side to begin with)
+  return
+    std::vector<std::pair<unsigned int, unsigned int> > ();
 }
 
 
@@ -673,14 +661,10 @@ template <int dim, int spacedim>
 FiniteElementDomination::Domination
 FE_DGQ<dim, spacedim>::compare_for_face_domination (const FiniteElement<dim, spacedim> &fe_other) const
 {
-  // check whether both are discontinuous
-  // elements, see the description of
-  // FiniteElementDomination::Domination
-  if (dynamic_cast<const FE_DGQ<dim, spacedim>*>(&fe_other) != 0)
-    return FiniteElementDomination::no_requirements;
-
-  Assert (false, ExcNotImplemented());
-  return FiniteElementDomination::neither_element_dominates;
+  // this is a discontinuous element, so by definition there will
+  // be no constraints wherever this element comes together with
+  // any other kind of element
+  return FiniteElementDomination::no_requirements;
 }
 
 
index 99c65476d8d359ad2e4ea31de71d460d0867679b..b384c1a3a101cb5816a5701804ecdacc33f3a1b9 100644 (file)
@@ -683,6 +683,17 @@ hp_vertex_dof_identities (const FiniteElement<dim,spacedim> &fe_other) const
       // equivalencies to be recorded
       return std::vector<std::pair<unsigned int, unsigned int> > ();
     }
+  else if (fe_other.dofs_per_face == 0)
+    {
+      // if the other element has no elements on faces at all,
+      // then it would be impossible to enforce any kind of
+      // continuity even if we knew exactly what kind of element
+      // we have -- simply because the other element declares
+      // that it is discontinuous because it has no DoFs on
+      // its faces. in that case, just state that we have no
+      // constraints to declare
+      return std::vector<std::pair<unsigned int, unsigned int> > ();
+    }
   else
     {
       Assert (false, ExcNotImplemented());
@@ -734,6 +745,17 @@ hp_line_dof_identities (const FiniteElement<dim,spacedim> &fe_other) const
       // equivalencies to be recorded
       return std::vector<std::pair<unsigned int, unsigned int> > ();
     }
+  else if (fe_other.dofs_per_face == 0)
+    {
+      // if the other element has no elements on faces at all,
+      // then it would be impossible to enforce any kind of
+      // continuity even if we knew exactly what kind of element
+      // we have -- simply because the other element declares
+      // that it is discontinuous because it has no DoFs on
+      // its faces. in that case, just state that we have no
+      // constraints to declare
+      return std::vector<std::pair<unsigned int, unsigned int> > ();
+    }
   else
     {
       Assert (false, ExcNotImplemented());
@@ -789,6 +811,17 @@ hp_quad_dof_identities (const FiniteElement<dim,spacedim>        &fe_other) cons
       // equivalencies to be recorded
       return std::vector<std::pair<unsigned int, unsigned int> > ();
     }
+  else if (fe_other.dofs_per_face == 0)
+    {
+      // if the other element has no elements on faces at all,
+      // then it would be impossible to enforce any kind of
+      // continuity even if we knew exactly what kind of element
+      // we have -- simply because the other element declares
+      // that it is discontinuous because it has no DoFs on
+      // its faces. in that case, just state that we have no
+      // constraints to declare
+      return std::vector<std::pair<unsigned int, unsigned int> > ();
+    }
   else
     {
       Assert (false, ExcNotImplemented());
diff --git a/tests/hp/q_vs_dgq.cc b/tests/hp/q_vs_dgq.cc
new file mode 100644 (file)
index 0000000..abfdf5a
--- /dev/null
@@ -0,0 +1,112 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2014 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.
+//
+// ---------------------------------------------------------------------
+
+
+
+// test by Korosh Taebi: check that we can gave Q(p) and DGQ(r) in the same
+// mesh
+
+
+#include "../tests.h"
+#include <deal.II/base/logstream.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/dofs/dof_accessor.h>
+#include <deal.II/hp/dof_handler.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_dgq.h>
+
+
+namespace Step27
+{
+  using namespace dealii;
+
+  template <int dim>
+  class MixedFECollection
+  {
+  public:
+    MixedFECollection ();
+    ~MixedFECollection ();
+
+    void run ();
+
+  private:
+
+    Triangulation<dim>          triangulation;
+    hp::DoFHandler<dim>      dof_handler;
+    hp::FECollection<dim>    fe_collection;
+
+  };
+
+  template <int dim>
+  MixedFECollection<dim>::MixedFECollection ()
+  :
+  dof_handler (triangulation)
+  {
+
+  }
+
+  template <int dim>
+  MixedFECollection<dim>::~MixedFECollection ()
+  {
+    dof_handler.clear ();
+  }
+
+  template <int dim>
+  void MixedFECollection<dim>::run ()
+  {
+    // add two a CG and a DG finite element object to fe_collection
+    fe_collection.push_back (FE_Q<dim>(1));
+    fe_collection.push_back (FE_DGQ<dim>(1));
+    deallog << " fe_collection size = " << fe_collection.size() << std::endl;
+
+    // produce a simple grid with 4 cells
+    GridGenerator::hyper_cube (triangulation, 0, 1);
+    triangulation.refine_global (1);
+
+    // looping over all cells and assigning the FE_DG object to the first cell
+    typename hp::DoFHandler<dim>::active_cell_iterator
+      cell = dof_handler.begin_active(),
+      endc = dof_handler.end();
+    for (unsigned int counter = 0; cell!=endc; ++cell, counter ++)
+      if (counter == 0)
+       {
+         cell->set_active_fe_index (cell->active_fe_index() + 1);
+       }
+
+    dof_handler.distribute_dofs (fe_collection);
+
+    deallog << "   Number of active cells:       " << triangulation.n_active_cells() << std::endl
+           << "   Number of degrees of freedom: " << dof_handler.n_dofs() << std::endl;
+  }
+}
+
+
+
+int main ()
+{
+  std::ofstream logfile("output");
+  logfile.precision(2);
+
+  deallog.attach(logfile);
+  deallog.depth_console(0);
+  deallog.threshold_double(1.e-10);
+
+  Step27::MixedFECollection<1>().run();
+  Step27::MixedFECollection<2>().run();
+  Step27::MixedFECollection<3>().run();
+
+  deallog << "OK" << std::endl;
+}
diff --git a/tests/hp/q_vs_dgq.output b/tests/hp/q_vs_dgq.output
new file mode 100644 (file)
index 0000000..16ded5a
--- /dev/null
@@ -0,0 +1,11 @@
+
+DEAL:: fe_collection size = 2
+DEAL::   Number of active cells:       2
+DEAL::   Number of degrees of freedom: 4
+DEAL:: fe_collection size = 2
+DEAL::   Number of active cells:       4
+DEAL::   Number of degrees of freedom: 12
+DEAL:: fe_collection size = 2
+DEAL::   Number of active cells:       8
+DEAL::   Number of degrees of freedom: 34
+DEAL::OK

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.