Implement function to check for MatrixFree support

diff --git a/doc/news/changes/minor/20170111DanielArndt b/doc/news/changes/minor/20170111DanielArndt
new file mode 100644 (file)
index 0000000..214782d
--- /dev/null
+++ b/doc/news/changes/minor/20170111DanielArndt
@@ -0,0 +1,3 @@
+New: MatrixFree allows to be asked if a certain FiniteElement is supported.
+<br>
+(Daniel Arndt, 2017/01/11)

diff --git a/include/deal.II/matrix_free/matrix_free.h b/include/deal.II/matrix_free/matrix_free.h
index a3e0d1d37bf1f728b98b4882cbcd2eeb633eabb4..0bf111a6cb795193d66a0691f4829b25b6c214bf 100644 (file)
--- a/include/deal.II/matrix_free/matrix_free.h
+++ b/include/deal.II/matrix_free/matrix_free.h
@@ -639,6 +639,13 @@ public:
    * @name 4: General information
    */
   //@{
+  /**
+   * Returns whether a given FiniteElement @p fe is supported by this class.
+   */
+  template <int spacedim>
+  static
+  bool is_supported (const FiniteElement<dim, spacedim> &fe);
+
   /**
    * Return the number of different DoFHandlers specified at initialization.
    */

diff --git a/include/deal.II/matrix_free/matrix_free.templates.h b/include/deal.II/matrix_free/matrix_free.templates.h
index 67e5a933440fabb044e03b4fdc4c3b128702f0d8..e30047393b833eb4ba233a2260cb1c326a43a8ef 100644 (file)
--- a/include/deal.II/matrix_free/matrix_free.templates.h
+++ b/include/deal.II/matrix_free/matrix_free.templates.h
@@ -331,6 +331,38 @@ internal_reinit(const Mapping<dim>                            &mapping,
 }
 
 
+template <int dim, typename Number>
+template <int spacedim>
+bool
+MatrixFree<dim, Number>::is_supported(const FiniteElement<dim, spacedim> &fe)
+{
+  if (dim != spacedim)
+    return false;
+
+  // first check for degree, number of base_elemnt and number of its components
+  if (fe.degree==0 || fe.n_base_elements()!=1)
+    return false;
+
+  const FiniteElement<dim, spacedim> *fe_ptr = &(fe.base_element(0));
+  if (fe_ptr->n_components() != 1)
+    return false;
+
+  // then check of the base element is supported
+  if (dynamic_cast<const FE_Poly<TensorProductPolynomials<dim>,dim,spacedim>*>(fe_ptr)!=0)
+    return true;
+  if (dynamic_cast<const FE_Poly<TensorProductPolynomials<dim,
+      Polynomials::PiecewisePolynomial<double> >,dim,spacedim>*>(fe_ptr)!=0)
+    return true;
+  if (dynamic_cast<const FE_DGP<dim, spacedim>*>(fe_ptr)!=0)
+    return true;
+  if (dynamic_cast<const FE_Q_DG0<dim, spacedim>*>(fe_ptr)!=0)
+    return true;
+
+  // if the base element is not in the above list it is not supported
+  return false;
+}
+
+
 
 namespace internal
 {

diff --git a/include/deal.II/numerics/vector_tools.templates.h b/include/deal.II/numerics/vector_tools.templates.h
index 920d959b7164afcfe3033a8704e0218dfdeb6fbb..0b64ae5dc32bf94aa75662895e9c7f56871f6f62 100644 (file)
--- a/include/deal.II/numerics/vector_tools.templates.h
+++ b/include/deal.II/numerics/vector_tools.templates.h
@@ -1113,34 +1113,15 @@ namespace VectorTools
                   const bool                    project_to_boundary_first)
     {
       // If we can, use the matrix-free implementation
-      bool use_matrix_free = true;
+      bool use_matrix_free = MatrixFree<dim, typename VectorType::value_type>::is_supported(dof.get_fe());
+
       // enforce_zero_boundary and project_to_boundary_first
-      // are not yet supported
-      if (enforce_zero_boundary || project_to_boundary_first)
+      // are not yet supported.
+      // We have explicit instantiations only if
+      // the number of components and the degree is not too high.
+      if (enforce_zero_boundary || project_to_boundary_first
+          || dof.get_fe().degree>8 || dof.get_fe().n_components()>4)
         use_matrix_free = false;
-      else
-        {
-          // Find out if the FiniteElement is supported
-          // This is based on matrix_free/shape_info.templates.h
-          // and matrix_free/matrix_free.templates.h
-          if (dof.get_fe().degree==0 || dof.get_fe().n_base_elements()!=1
-              || dof.get_fe().degree>8 || dof.get_fe().n_components()>4)
-            use_matrix_free = false;
-          else
-            {
-              const FiniteElement<dim> *fe_ptr = &dof.get_fe().base_element(0);
-              if (fe_ptr->n_components() != 1)
-                use_matrix_free = false;
-              else if ((dynamic_cast<const FE_Poly<TensorProductPolynomials<dim>,dim,dim>*>(fe_ptr)==0)
-                       && (dynamic_cast<const FE_Poly<TensorProductPolynomials<dim,
-                           Polynomials::PiecewisePolynomial<double> >,dim,dim>*>(fe_ptr)==0)
-                       &&(dynamic_cast<const FE_DGP<dim>*>(fe_ptr)==0)
-                       &&(dynamic_cast<const FE_Q_DG0<dim>*>(fe_ptr)==0)
-                       &&(dynamic_cast<const FE_Q<dim>*>(fe_ptr)==0)
-                       &&(dynamic_cast<const FE_DGQ<dim>*>(fe_ptr)==0))
-                use_matrix_free = false;
-            }
-        }
 
       if (use_matrix_free)
         project_matrix_free_component (mapping, dof, constraints, quadrature,

diff --git a/source/matrix_free/matrix_free.inst.in b/source/matrix_free/matrix_free.inst.in
index aaa35c70c8dc9328926dd1d51855772229761f77..fe24f54f217f58619fc0c38f1637feae75695f13 100644 (file)
--- a/source/matrix_free/matrix_free.inst.in
+++ b/source/matrix_free/matrix_free.inst.in
@@ -1,6 +1,6 @@
 // ---------------------------------------------------------------------
 //
-// Copyright (C) 2010 - 2016 by the deal.II authors
+// Copyright (C) 2010 - 2017 by the deal.II authors
 //
 // This file is part of the deal.II library.
 //
@@ -41,3 +41,14 @@ for (deal_II_dimension : DIMENSIONS)
                         <deal_II_dimension,deal_II_dimension> &, const unsigned int);
 #endif
 }
+
+for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension :  SPACE_DIMENSIONS)
+{
+#if deal_II_dimension <= deal_II_space_dimension
+    template bool MatrixFree<deal_II_dimension, double>::
+    is_supported(const FiniteElement<deal_II_dimension, deal_II_space_dimension> &);
+
+    template bool MatrixFree<deal_II_dimension, float>::
+    is_supported(const FiniteElement<deal_II_dimension, deal_II_space_dimension> &);
+#endif
+}

diff --git a/tests/matrix_free/is_supported_01.cc b/tests/matrix_free/is_supported_01.cc
new file mode 100644 (file)
index 0000000..19a368a
--- /dev/null
+++ b/tests/matrix_free/is_supported_01.cc
@@ -0,0 +1,90 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2017 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 the ouput of MatrixFree::is_supported for various FiniteElements
+
+#include "../tests.h"
+#include <fstream>
+#include <deal.II/base/logstream.h>
+
+#include <deal.II/fe/fe_q_iso_q1.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_q_dg0.h>
+#include <deal.II/fe/fe_q_bubbles.h>
+#include <deal.II/fe/fe_bernstein.h>
+#include <deal.II/fe/fe_q_hierarchical.h>
+#include <deal.II/fe/fe_dgq.h>
+#include <deal.II/fe/fe_dgp.h>
+#include <deal.II/fe/fe_dgp_monomial.h>
+#include <deal.II/fe/fe_dgp_nonparametric.h>
+#include <deal.II/fe/fe_dg_vector.h>
+#include <deal.II/fe/fe_nedelec.h>
+#include <deal.II/fe/fe_abf.h>
+#include <deal.II/fe/fe_bdm.h>
+#include <deal.II/fe/fe_raviart_thomas.h>
+#include <deal.II/fe/fe_rannacher_turek.h>
+#include <deal.II/fe/fe_nothing.h>
+#include <deal.II/fe/fe_face.h>
+#include <deal.II/fe/fe_p1nc.h>
+#include <deal.II/fe/fe_trace.h>
+
+#include <deal.II/fe/fe_tools.h>
+#include <deal.II/matrix_free/matrix_free.h>
+
+
+template <int dim, int spacedim>
+void print(const FiniteElement<dim, spacedim> &fe)
+{
+  deallog << fe.get_name()  << " supported by MatrixFree: "
+          << std::boolalpha << MatrixFree<dim>::is_supported(fe)
+          << std::endl;
+}
+
+
+int main ()
+{
+  initlog();
+
+  print<2,2> (FE_ABF<2>(0));
+  print<2,2> (FE_BDM<2>(1));
+  print<2,2> (FE_Bernstein<2>(1));
+  print<2,2> (FE_DGBDM<2>(1));
+  print<2,2> (FE_DGNedelec<2>(1));
+  print<1,2> (FE_DGP<1,2>(1));
+  print<2,2> (FE_DGP<2>(1));
+  print<2,2> (FE_DGPMonomial<2>(1));
+  print<2,2> (FE_DGPNonparametric<2>(1));
+  print<2,2> (FE_DGRaviartThomas<2>(1));
+  print<1,2> (FE_DGQ<1,2>(1));
+  print<2,2> (FE_DGQ<2>(1));
+  print<2,2> (FE_DGQArbitraryNodes<2>(QGauss<1>(2)));
+  print<2,2> (FE_FaceP<2>(1));
+  print<2,2> (FE_FaceQ<2>(1));
+  print<2,2> (FE_P1NC());
+  print<2,2> (FE_Nedelec<2>(1));
+  print<2,2> (FE_Nothing<2>());
+  print<1,2> (FE_Q<1,2>(1));
+  print<2,2> (FE_Q<2>(1));
+  print<2,2> (FE_Q_Bubbles<2>(1));
+  print<2,2> (FE_Q_DG0<2>(1));
+  print<2,2> (FE_Q_Hierarchical<2>(1));
+  print<2,2> (FE_Q_iso_Q1<2>(1));
+  print<2,2> (FE_RannacherTurek<2>(0));
+  print<2,2> (FE_RaviartThomas<2>(1));
+  print<2,2> (FE_RaviartThomasNodal<2>(1));
+  print<2,2> (FE_TraceQ<2>(1));
+
+  return 0;
+}

diff --git a/tests/matrix_free/is_supported_01.output b/tests/matrix_free/is_supported_01.output
new file mode 100644 (file)
index 0000000..2601ac3
--- /dev/null
+++ b/tests/matrix_free/is_supported_01.output
@@ -0,0 +1,29 @@
+
+DEAL::FE_ABF<2>(0) supported by MatrixFree: false
+DEAL::FE_BDM<2>(1) supported by MatrixFree: false
+DEAL::FE_Bernstein<2>(1) supported by MatrixFree: true
+DEAL::FE_DGBDM<2>(1) supported by MatrixFree: false
+DEAL::FE_DGNedelec<2>(1) supported by MatrixFree: false
+DEAL::FE_DGP<1,2>(1) supported by MatrixFree: false
+DEAL::FE_DGP<2>(1) supported by MatrixFree: true
+DEAL::FE_DGPMonomial<2>(1) supported by MatrixFree: false
+DEAL::FE_DGPNonparametric<2>(1) supported by MatrixFree: false
+DEAL::FE_DGRaviartThomas<2>(1) supported by MatrixFree: false
+DEAL::FE_DGQ<1,2>(1) supported by MatrixFree: false
+DEAL::FE_DGQ<2>(1) supported by MatrixFree: true
+DEAL::FE_DGQArbitraryNodes<2>(QGauss(2)) supported by MatrixFree: true
+DEAL::FE_FaceP<2>(1) supported by MatrixFree: false
+DEAL::FE_FaceQ<2>(1) supported by MatrixFree: false
+DEAL::FE_P1NC supported by MatrixFree: false
+DEAL::FE_Nedelec<2>(1) supported by MatrixFree: false
+DEAL::FE_Nothing<2>() supported by MatrixFree: false
+DEAL::FE_Q<1,2>(1) supported by MatrixFree: false
+DEAL::FE_Q<2>(1) supported by MatrixFree: true
+DEAL::FE_Q_Bubbles<2>(1) supported by MatrixFree: false
+DEAL::FE_Q_DG0<2>(1) supported by MatrixFree: true
+DEAL::FE_Q_Hierarchical<2>(1) supported by MatrixFree: true
+DEAL::FE_Q_iso_Q1<2>(1) supported by MatrixFree: true
+DEAL::FE_RannacherTurek<2>(0, 2) supported by MatrixFree: false
+DEAL::FE_RaviartThomas<2>(1) supported by MatrixFree: false
+DEAL::FE_RaviartThomasNodal<2>(1) supported by MatrixFree: false
+DEAL::FE_TraceQ<2>(1) supported by MatrixFree: false