]> https://gitweb.dealii.org/ - dealii.git/commitdiff
extend VectorTools::project() to take fe_component 5874/head
authorDenis Davydov <davydden@gmail.com>
Wed, 7 Feb 2018 13:24:37 +0000 (14:24 +0100)
committerDenis Davydov <davydden@gmail.com>
Thu, 8 Feb 2018 08:28:01 +0000 (09:28 +0100)
doc/news/changes/minor/20180207DenisDavydov [new file with mode: 0644]
include/deal.II/numerics/vector_tools.h
include/deal.II/numerics/vector_tools.templates.h
source/numerics/vector_tools_project_qpmf.inst.in
tests/numerics/project_parallel_qpmf_common.h
tests/numerics/project_parallel_qpmf_q_03.cc [new file with mode: 0644]
tests/numerics/project_parallel_qpmf_q_03.with_p4est=true.mpirun=3.output [new file with mode: 0644]

diff --git a/doc/news/changes/minor/20180207DenisDavydov b/doc/news/changes/minor/20180207DenisDavydov
new file mode 100644 (file)
index 0000000..4c1769f
--- /dev/null
@@ -0,0 +1,4 @@
+Extend VectorTools::project() function with MatrixFree quadrature data to
+optionally take the finite element component index.
+<br>
+(Denis Davydov, 2018/02/07)
index 55eb3a43316aa2308c85900e5e344b376a4af869..d0c121b2cde7908ce7737c4500f317e13d00fbac 100644 (file)
@@ -919,13 +919,18 @@ namespace VectorTools
    * @endcode
    * where <code>qp_data</code> is a an object of type Table<2, VectorizedArray<double> >,
    * which stores quadrature point data.
+   *
+   * @p fe_component allow to additionally specify which component of @p data
+   * to use in case it was constructed with an <code>std::vector<const DoFHandler<dim>*></code>
+   * in will be internally used in constructor of FEEvaluation object.
    */
   template <int dim, typename VectorType>
   void project (std::shared_ptr<const MatrixFree<dim,typename VectorType::value_type> > data,
                 const ConstraintMatrix &constraints,
                 const unsigned int n_q_points_1d,
                 const std::function< VectorizedArray<typename VectorType::value_type> (const unsigned int, const unsigned int)> &func,
-                VectorType &vec_result);
+                VectorType &vec_result,
+                const unsigned int fe_component = 0);
 
   /**
    * Same as above but for <code>n_q_points_1d = matrix_free.get_dof_handler().get_fe().degree+1</code>.
@@ -934,7 +939,8 @@ namespace VectorTools
   void project (std::shared_ptr<const MatrixFree<dim,typename VectorType::value_type> > data,
                 const ConstraintMatrix &constraints,
                 const std::function< VectorizedArray<typename VectorType::value_type> (const unsigned int, const unsigned int)> &func,
-                VectorType &vec_result);
+                VectorType &vec_result,
+                const unsigned int fe_component = 0);
 
   /**
    * Compute Dirichlet boundary conditions.  This function makes up a map of
index 1e68a740e08b9c5022da3433bb8810ed478ab932..d07d0bbdce3a0af24296840dfc36919db4857e27 100644 (file)
@@ -1386,9 +1386,10 @@ namespace VectorTools
     void project_parallel (std::shared_ptr<const MatrixFree<dim,typename VectorType::value_type> > matrix_free,
                            const ConstraintMatrix &constraints,
                            const std::function< VectorizedArray<typename VectorType::value_type> (const unsigned int, const unsigned int)> &func,
-                           VectorType &vec_result)
+                           VectorType &vec_result,
+                           const unsigned int fe_component)
     {
-      const DoFHandler<dim,spacedim> &dof = matrix_free->get_dof_handler();
+      const DoFHandler<dim,spacedim> &dof = matrix_free->get_dof_handler(fe_component);
 
       typedef typename VectorType::value_type Number;
       Assert (dof.get_fe(0).n_components() == 1,
@@ -1402,20 +1403,20 @@ namespace VectorTools
 
       typedef MatrixFreeOperators::MassOperator<dim, fe_degree, n_q_points_1d, 1, LinearAlgebra::distributed::Vector<Number> > MatrixType;
       MatrixType mass_matrix;
-      mass_matrix.initialize(matrix_free);
+      mass_matrix.initialize(matrix_free, {{fe_component}});
       mass_matrix.compute_diagonal();
 
       typedef LinearAlgebra::distributed::Vector<Number> LocalVectorType;
       LocalVectorType vec, rhs, inhomogeneities;
-      matrix_free->initialize_dof_vector(vec);
-      matrix_free->initialize_dof_vector(rhs);
-      matrix_free->initialize_dof_vector(inhomogeneities);
+      matrix_free->initialize_dof_vector(vec,fe_component);
+      matrix_free->initialize_dof_vector(rhs,fe_component);
+      matrix_free->initialize_dof_vector(inhomogeneities,fe_component);
       constraints.distribute(inhomogeneities);
       inhomogeneities*=-1.;
 
       // assemble right hand side:
       {
-        FEEvaluation<dim,fe_degree,n_q_points_1d,1,Number> fe_eval(*matrix_free);
+        FEEvaluation<dim,fe_degree,n_q_points_1d,1,Number> fe_eval(*matrix_free, fe_component);
         const unsigned int n_cells = matrix_free->n_macro_cells();
         const unsigned int n_q_points = fe_eval.n_q_points;
 
@@ -1489,28 +1490,28 @@ namespace VectorTools
                 const ConstraintMatrix &constraints,
                 const unsigned int n_q_points_1d,
                 const std::function< VectorizedArray<typename VectorType::value_type> (const unsigned int, const unsigned int)> &func,
-                VectorType &vec_result)
+                VectorType &vec_result,
+                const unsigned int fe_component)
   {
-    (void) n_q_points_1d;
-    const unsigned int fe_degree = matrix_free->get_dof_handler().get_fe().degree;
-
-    Assert (fe_degree+1 == n_q_points_1d,
-            ExcNotImplemented());
+    const unsigned int fe_degree = matrix_free->get_dof_handler(fe_component).get_fe().degree;
 
-    switch (fe_degree)
-      {
-      case 1:
-        project_parallel<dim,VectorType,dim,1,2> (matrix_free,constraints,func,vec_result);
-        break;
-      case 2:
-        project_parallel<dim,VectorType,dim,2,3> (matrix_free,constraints,func,vec_result);
-        break;
-      case 3:
-        project_parallel<dim,VectorType,dim,3,4> (matrix_free,constraints,func,vec_result);
-        break;
-      default:
-        project_parallel<dim,VectorType,dim,-1,0> (matrix_free,constraints,func,vec_result);
-      }
+    if (fe_degree+1 == n_q_points_1d)
+      switch (fe_degree)
+        {
+        case 1:
+          project_parallel<dim,VectorType,dim,1,2> (matrix_free,constraints,func,vec_result,fe_component);
+          break;
+        case 2:
+          project_parallel<dim,VectorType,dim,2,3> (matrix_free,constraints,func,vec_result,fe_component);
+          break;
+        case 3:
+          project_parallel<dim,VectorType,dim,3,4> (matrix_free,constraints,func,vec_result,fe_component);
+          break;
+        default:
+          project_parallel<dim,VectorType,dim,-1,0> (matrix_free,constraints,func,vec_result,fe_component);
+        }
+    else
+      project_parallel<dim,VectorType,dim,-1,0> (matrix_free,constraints,func,vec_result,fe_component);
   }
 
 
@@ -1519,13 +1520,15 @@ namespace VectorTools
   void project (std::shared_ptr<const MatrixFree<dim,typename VectorType::value_type> > matrix_free,
                 const ConstraintMatrix &constraints,
                 const std::function< VectorizedArray<typename VectorType::value_type> (const unsigned int, const unsigned int)> &func,
-                VectorType &vec_result)
+                VectorType &vec_result,
+                const unsigned int fe_component)
   {
     project (matrix_free,
              constraints,
-             matrix_free->get_dof_handler().get_fe().degree+1,
+             matrix_free->get_dof_handler(fe_component).get_fe().degree+1,
              func,
-             vec_result);
+             vec_result,
+             fe_component);
   }
 
 
index 3b9c0968c6af48505dfd7ae7f096003ee6036297..4db6e507cbd70670555ea2d6aba93d5905cbe9f1 100644 (file)
@@ -24,7 +24,8 @@ for (VEC: REAL_NONBLOCK_VECTORS; deal_II_dimension : DIMENSIONS)
         std::shared_ptr<const MatrixFree<deal_II_dimension,VEC::value_type> > matrix_free,
         const ConstraintMatrix &constraints,
         const std::function< VectorizedArray<VEC::value_type> (const unsigned int, const unsigned int)> &,
-        VEC &);
+        VEC &,
+        const unsigned int);
 
     template
     void project<deal_II_dimension, VEC>(
@@ -32,7 +33,8 @@ for (VEC: REAL_NONBLOCK_VECTORS; deal_II_dimension : DIMENSIONS)
         const ConstraintMatrix &constraints,
         const unsigned int,
         const std::function< VectorizedArray<VEC::value_type> (const unsigned int, const unsigned int)> &,
-        VEC &);
+        VEC &,
+        const unsigned int);
 
     \}
 
index 79c1ec6143bb24d86887433189c6912931b88b9a..d83496e977e51e970d64afe6fac54a3406a5e19e 100644 (file)
@@ -48,6 +48,9 @@
 #include <fstream>
 #include <vector>
 
+/**
+ * A function of polynomial degree @p q with n_components.
+ */
 template <int dim>
 class F :  public Function<dim>
 {
@@ -90,27 +93,44 @@ private:
 };
 
 
+/**
+ * A general function to perform projection of polynomial function of degrees [0,p]
+ * onto the FE space given by @p fe and @p triangulation .
+ */
 template <int fe_degree, int n_q_points_1d, int dim>
 void do_project (const parallel::distributed::Triangulation<dim> &triangulation,
-                 const FiniteElement<dim> &fe,
-                 const unsigned int        p)
+                 const std::vector<const FiniteElement<dim>*> &fes,
+                 const unsigned int        p,
+                 const unsigned int        fe_index = 0)
 {
-  AssertThrow(fe.n_components() ==1,
+  // we only project scalar value function
+  AssertThrow(fes[fe_index]->n_components() ==1,
               ExcNotImplemented());
-  DoFHandler<dim>        dof_handler(triangulation);
-  dof_handler.distribute_dofs (fe);
 
   deallog << "n_cells=" << triangulation.n_global_active_cells() << std::endl;
-  deallog << "n_dofs=" << dof_handler.n_dofs() << std::endl;
 
-  IndexSet locally_relevant_dofs;
-  DoFTools::extract_locally_relevant_dofs (dof_handler,
-                                           locally_relevant_dofs);
+  std::vector<std::shared_ptr<DoFHandler<dim>>> dof_handlers(fes.size());
+  std::vector<IndexSet> locally_relevant_dofs(fes.size());
+  std::vector<ConstraintMatrix> constraints(fes.size());
 
-  ConstraintMatrix constraints;
-  constraints.reinit (locally_relevant_dofs);
-  DoFTools::make_hanging_node_constraints (dof_handler, constraints);
-  constraints.close ();
+  std::vector<const DoFHandler<dim>*> dof_handlers_mf(fes.size());
+  std::vector<const ConstraintMatrix *> constraints_mf(fes.size());
+  for (unsigned int i = 0; i < fes.size(); ++i)
+    {
+      dof_handlers[i] = std::make_shared<DoFHandler<dim>>(triangulation);
+      dof_handlers[i]->distribute_dofs (*fes[i]);
+      deallog << "n_dofs=" << dof_handlers[i]->n_dofs() << std::endl;
+
+      DoFTools::extract_locally_relevant_dofs (*dof_handlers[i],
+                                               locally_relevant_dofs[i]);
+
+      constraints[i].reinit (locally_relevant_dofs[i]);
+      DoFTools::make_hanging_node_constraints (*dof_handlers[i], constraints[i]);
+      constraints[i].close ();
+
+      dof_handlers_mf[i] = dof_handlers[i].get();
+      constraints_mf[i]  = &constraints[i];
+    }
 
   QGauss<1> quadrature_formula_1d(n_q_points_1d);
 
@@ -120,16 +140,16 @@ void do_project (const parallel::distributed::Triangulation<dim> &triangulation,
   additional_data.tasks_parallel_scheme = MatrixFree<dim,double>::AdditionalData::partition_color;
   additional_data.mapping_update_flags = update_values | update_JxW_values | update_quadrature_points;
   std::shared_ptr<MatrixFree<dim,double> >  data(new MatrixFree<dim,double> ());
-  data->reinit (dof_handler, constraints, quadrature_formula_1d, additional_data);
+  data->reinit (dof_handlers_mf, constraints_mf, quadrature_formula_1d, additional_data);
 
   for (unsigned int q=0; q<=p; ++q)
     {
       // setup quadrature data:
-      F<dim> function(q, fe.n_components());
+      F<dim> function(q, fes[fe_index]->n_components());
 
       // initialize a quadrature data
       {
-        FEEvaluation<dim,fe_degree,n_q_points_1d,1,double> fe_eval(*data);
+        FEEvaluation<dim,fe_degree,n_q_points_1d,1,double> fe_eval(*data,fe_index);
         const unsigned int n_cells = data->n_macro_cells();
         const unsigned int n_q_points = fe_eval.n_q_points;
 
@@ -143,13 +163,14 @@ void do_project (const parallel::distributed::Triangulation<dim> &triangulation,
       }
 
       LinearAlgebra::distributed::Vector<double> field;
-      data->initialize_dof_vector(field);
+      data->initialize_dof_vector(field,fe_index);
       VectorTools::project<dim,LinearAlgebra::distributed::Vector<double> >
       (data,
-       constraints,
+       constraints[fe_index],
        n_q_points_1d,
        [=] (const unsigned int cell, const unsigned int q) -> VectorizedArray<double> { return qp_data(cell,q); },
-       field);
+       field,
+       fe_index);
 
       field.update_ghost_values();
 
@@ -159,16 +180,16 @@ void do_project (const parallel::distributed::Triangulation<dim> &triangulation,
       double L2_norm = 0.;
       {
         QGauss<dim> quadrature_formula_error(std::max(p,q)+1);
-        FEValues<dim> fe_values (fe, quadrature_formula_error,
+        FEValues<dim> fe_values (*fes[fe_index], quadrature_formula_error,
                                  update_values  | update_quadrature_points | update_JxW_values);
 
-        const unsigned int   dofs_per_cell = fe.dofs_per_cell;
+        const unsigned int   dofs_per_cell = fes[fe_index]->dofs_per_cell;
         const unsigned int   n_q_points    = quadrature_formula_error.size();
         std::vector<double>  values (n_q_points);
 
         typename DoFHandler<dim>::active_cell_iterator
-        cell = dof_handler.begin_active(),
-        endc = dof_handler.end();
+        cell = dof_handlers[fe_index]->begin_active(),
+        endc = dof_handlers[fe_index]->end();
         for (; cell!=endc; ++cell)
           if (cell->is_locally_owned())
             {
@@ -186,7 +207,7 @@ void do_project (const parallel::distributed::Triangulation<dim> &triangulation,
       L2_norm = Utilities::MPI::sum(L2_norm, MPI_COMM_WORLD);
       L2_norm = std::sqrt(L2_norm);
 
-      deallog << fe.get_name() << ", P_" << q
+      deallog << fes[fe_index]->get_name() << ", P_" << q
               << ", rel. error=" << L2_norm / field_l2_norm
               << std::endl;
 
@@ -212,7 +233,7 @@ void test_no_hanging_nodes (const FiniteElement<dim> &fe,
   GridGenerator::hyper_cube (triangulation);
   triangulation.refine_global (3);
 
-  do_project<fe_degree, n_q_points_1d, dim> (triangulation, fe, p);
+  do_project<fe_degree, n_q_points_1d, dim> (triangulation, {{&fe}}, p);
 }
 
 
@@ -230,6 +251,23 @@ void test_with_hanging_nodes (const FiniteElement<dim> &fe,
   triangulation.execute_coarsening_and_refinement ();
   triangulation.refine_global (1);
 
-  do_project<fe_degree, n_q_points_1d, dim> (triangulation, fe, p);
+  do_project<fe_degree, n_q_points_1d, dim> (triangulation, {{&fe}}, p);
 }
 
+
+// same as above but for multiple fes
+template <int fe_degree, int n_q_points_1d, int dim>
+void test_with_hanging_nodes (const std::vector<const FiniteElement<dim>*> &fes,
+                              const unsigned int        p,
+                              const unsigned int fe_index)
+{
+  parallel::distributed::Triangulation<dim> triangulation(MPI_COMM_WORLD);
+  GridGenerator::hyper_cube (triangulation);
+  triangulation.refine_global (2);
+  if (triangulation.begin_active()->is_locally_owned())
+    triangulation.begin_active()->set_refine_flag();
+  triangulation.execute_coarsening_and_refinement ();
+  triangulation.refine_global (1);
+
+  do_project<fe_degree, n_q_points_1d, dim> (triangulation, fes, p, fe_index);
+}
diff --git a/tests/numerics/project_parallel_qpmf_q_03.cc b/tests/numerics/project_parallel_qpmf_q_03.cc
new file mode 100644 (file)
index 0000000..fadf832
--- /dev/null
@@ -0,0 +1,43 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2018 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.
+//
+// ---------------------------------------------------------------------
+
+
+// check VectorTools::project_parallel() for matrix-free quadrature data for
+// FE_Q on a mesh with hanging nodes and specified fe_component
+
+#include "../tests.h"
+#include "project_parallel_qpmf_common.h"
+
+template <int dim>
+void test ()
+{
+  FE_Q<dim> fe1(1);
+  FE_Q<dim> fe2(2);
+  FE_Q<dim> fe4(4);
+  test_with_hanging_nodes<1, 2,dim> ({{&fe1,&fe2}}, 1, 0); // take first
+  test_with_hanging_nodes<2, 3,dim> ({{&fe1,&fe2}}, 2, 1); // take second
+  test_with_hanging_nodes<2, 6,dim> ({{&fe2,&fe4}}, 2, 0); // take first with non-standard n_q_points_1d
+}
+
+int main(int argc, char *argv[])
+{
+  Utilities::MPI::MPI_InitFinalize mpi_initialization (argc, argv,
+                                                       numbers::invalid_unsigned_int);
+
+  initlog();
+
+  test<2>();
+  test<3>();
+}
diff --git a/tests/numerics/project_parallel_qpmf_q_03.with_p4est=true.mpirun=3.output b/tests/numerics/project_parallel_qpmf_q_03.with_p4est=true.mpirun=3.output
new file mode 100644 (file)
index 0000000..9fd643d
--- /dev/null
@@ -0,0 +1,35 @@
+
+DEAL::n_cells=76
+DEAL::n_dofs=97
+DEAL::n_dofs=349
+DEAL::FE_Q<2>(1), P_0, rel. error=0
+DEAL::FE_Q<2>(1), P_1, rel. error=0
+DEAL::n_cells=76
+DEAL::n_dofs=97
+DEAL::n_dofs=349
+DEAL::FE_Q<2>(2), P_0, rel. error=0
+DEAL::FE_Q<2>(2), P_1, rel. error=0
+DEAL::FE_Q<2>(2), P_2, rel. error=0
+DEAL::n_cells=76
+DEAL::n_dofs=349
+DEAL::n_dofs=1309
+DEAL::FE_Q<2>(2), P_0, rel. error=0
+DEAL::FE_Q<2>(2), P_1, rel. error=0
+DEAL::FE_Q<2>(2), P_2, rel. error=0
+DEAL::n_cells=568
+DEAL::n_dofs=827
+DEAL::n_dofs=5559
+DEAL::FE_Q<3>(1), P_0, rel. error=0
+DEAL::FE_Q<3>(1), P_1, rel. error=0
+DEAL::n_cells=568
+DEAL::n_dofs=827
+DEAL::n_dofs=5559
+DEAL::FE_Q<3>(2), P_0, rel. error=0
+DEAL::FE_Q<3>(2), P_1, rel. error=0
+DEAL::FE_Q<3>(2), P_2, rel. error=0
+DEAL::n_cells=568
+DEAL::n_dofs=5559
+DEAL::n_dofs=40319
+DEAL::FE_Q<3>(2), P_0, rel. error=0
+DEAL::FE_Q<3>(2), P_1, rel. error=0
+DEAL::FE_Q<3>(2), P_2, rel. error=0

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.