]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Using the SolutionTransfer class with hp::DoFHandler and on meshes where some cells...
authorWolfgang Bangerth <bangerth@math.tamu.edu>
Fri, 29 Jun 2012 08:25:48 +0000 (08:25 +0000)
committerWolfgang Bangerth <bangerth@math.tamu.edu>
Fri, 29 Jun 2012 08:25:48 +0000 (08:25 +0000)
git-svn-id: https://svn.dealii.org/trunk@25660 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/doc/news/changes.h
deal.II/source/numerics/solution_transfer.cc
tests/hp/solution_transfer_02.cc [new file with mode: 0644]
tests/hp/solution_transfer_02/cmp/generic [new file with mode: 0644]

index b30926173351db32b58bbd0ec9c40efd2e8c105b..49d651cd180413b52742b509c0d0a8038054c94b 100644 (file)
@@ -47,6 +47,12 @@ used for boundary indicators.
 
 
 <ol>
+<li> Fixed: Using the SolutionTransfer class with hp::DoFHandler
+and on meshes where some cells are associated with a FE_Nothing element
+could result in an error. This is now fixed.
+<br>
+(Wolfgang Bangerth, 2012/06/29)
+
 <li> Fixed: The MappingQ1::transform_real_to_unit_cell function as
 well as the equivalent ones in derived classes sometimes get into
 trouble if they are asked to compute the preimage of this point
index cb6f005551e07dd8be21871df9a1f0fc61937899..f7799903067eecea4999d101cef9769145f2b5e1 100644 (file)
@@ -154,6 +154,19 @@ SolutionTransfer<dim, VECTOR, DH>::refine_interpolate(const VECTOR &in,
 
 namespace internal
 {
+                                  /**
+                                   * Generate a table that contains
+                                   * interpolation matrices between
+                                   * each combination of finite
+                                   * elements used in a DoFHandler of
+                                   * some kind. Since not all
+                                   * elements can be interpolated
+                                   * onto each other, the table may
+                                   * contain empty matrices for those
+                                   * combinations of elements for
+                                   * which no such interpolation is
+                                   * implemented.
+                                   */
   template <typename DH>
   void extract_interpolation_matrices (const DH&,
                                        Table<2,FullMatrix<double> > &)
@@ -170,7 +183,33 @@ namespace internal
         if (i != j)
           {
             matrices(i,j).reinit (fe[i].dofs_per_cell, fe[j].dofs_per_cell);
-            fe[i].get_interpolation_matrix (fe[j], matrices(i,j));
+
+                                            // see if we can get the
+                                            // interpolation matrices
+                                            // for this combination
+                                            // of elements. if not,
+                                            // reset the matrix sizes
+                                            // to zero to indicate
+                                            // that this particular
+                                            // combination isn't
+                                            // supported. this isn't
+                                            // an outright error
+                                            // right away since we
+                                            // may never need to
+                                            // actually interpolate
+                                            // between these two
+                                            // elements on actual
+                                            // cells; we simply have
+                                            // to trigger an error if
+                                            // someone actually tries
+           try
+             {
+               fe[i].get_interpolation_matrix (fe[j], matrices(i,j));
+             }
+           catch (const typename FiniteElement<dim,spacedim>::ExcInterpolationNotImplemented &)
+             {
+               matrices(i,j).reinit (0,0);
+             }
           }
   }
 
@@ -347,7 +386,50 @@ prepare_for_coarsening_and_refinement(const std::vector<VECTOR> &all_in)
                         {
                           tmp2.reinit (cell->child(most_general_child)->get_fe().dofs_per_cell,
                                        true);
-                          interpolation_hp (fe_ind_general, child_ind).vmult (tmp2, tmp);
+
+                                                          // if the
+                                                          // matrix
+                                                          // has size
+                                                          // 0 and
+                                                          // the
+                                                          // corresponding
+                                                          // elements
+                                                          // have
+                                                          // more
+                                                          // than
+                                                          // zero
+                                                          // DoFs,
+                                                          // then
+                                                          // this
+                                                          // means
+                                                          // that the
+                                                          // internal::extract_interpolation_matrices
+                                                          // function
+                                                          // above
+                                                          // couldn't
+                                                          // get the
+                                                          // interpolation
+                                                          // matrix
+                                                          // between
+                                                          // this
+                                                          // pair of
+                                                          // elements. since
+                                                          // we need
+                                                          // it here,
+                                                          // this is
+                                                          // an error
+                         if ((interpolation_hp(fe_ind_general, child_ind).m() == 0)
+                             &&
+                             (interpolation_hp(fe_ind_general, child_ind).n() == 0)
+                             &&
+                             ((tmp2.size() > 0) || (tmp.size() > 0)))
+                           Assert (false,
+                                   ExcMessage (std::string("No interpolation from element ")
+                                               + cell->child(child)->get_fe().get_name()
+                                               + " to element "
+                                               + cell->child(most_general_child)->get_fe().get_name()
+                                               + " is available, but it is needed here."));
+                          interpolation_hp(fe_ind_general, child_ind).vmult (tmp2, tmp);
                         }
                       else
                         tmp2.swap (tmp);
diff --git a/tests/hp/solution_transfer_02.cc b/tests/hp/solution_transfer_02.cc
new file mode 100644 (file)
index 0000000..33de200
--- /dev/null
@@ -0,0 +1,108 @@
+//----------------------------  solution_transfer_02.cc  ---------------------------
+//    fe_data_test.cc,v 1.14 2003/11/28 11:52:35 guido Exp
+//    Version:
+//
+//    Copyright (C) 2011, 2012 by the deal.II authors
+//
+//    This file is subject to QPL and may not be  distributed
+//    without copyright and license information. Please refer
+//    to the file deal.II/doc/license.html for the  text  and
+//    further information on this license.
+//
+//----------------------------  solution_transfer_02.cc  ---------------------------
+
+
+// SolutionTransfer wanted to compute interpolation matrices between
+// all pairs of elements used on a mesh in the hp case. unfortunately,
+// not all pairs are actually supported, e.g. between FE_Nothing and
+// FE_Q, but that shouldn't matter as long as these combinations are
+// never exercised on actual cells
+
+#include "../tests.h"
+#include <deal.II/base/function.h>
+#include <deal.II/base/logstream.h>
+#include <deal.II/lac/vector.h>
+
+#include <deal.II/grid/tria.h>
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_refinement.h>
+#include <deal.II/grid/tria_accessor.h>
+#include <deal.II/grid/tria_iterator.h>
+#include <deal.II/dofs/dof_accessor.h>
+#include <deal.II/dofs/dof_tools.h>
+#include <deal.II/numerics/solution_transfer.h>
+#include <deal.II/fe/fe_nothing.h>
+#include <deal.II/fe/fe_q.h>
+#include <fstream>
+#include <iostream>
+#include <vector>
+
+
+template <int dim>
+void transfer(std::ostream &out)
+{
+  Triangulation<dim> tria;
+  GridGenerator::hyper_cube(tria);
+  tria.refine_global(1);
+
+  hp::FECollection<dim> fe;
+  fe.push_back (FE_Q<dim> (1));
+  fe.push_back (FE_Nothing<dim>());
+
+                                  // create a DoFHandler on which we
+                                  // have both cells with FE_Q as
+                                  // well as FE_Nothing
+  hp::DoFHandler<dim> dof_handler(tria);
+  dof_handler.begin(0)->child(0)->set_active_fe_index(1);
+
+  Vector<double> solution;
+  ConstraintMatrix cm;
+  cm.close();
+
+  dof_handler.distribute_dofs (fe);
+  solution.reinit(dof_handler.n_dofs());
+
+  for (unsigned int i=0; i<solution.size(); ++i)
+    solution(i) = i;
+
+  SolutionTransfer<dim,Vector<double>,hp::DoFHandler<dim> > soltrans(dof_handler);
+
+  typename Triangulation<dim>::active_cell_iterator cell=tria.begin_active(),
+                                                   endc=tria.end();
+  ++cell;
+  ++cell;
+  for (; cell!=endc; ++cell)
+    cell->set_refine_flag();
+
+  Vector<double> old_solution=solution;
+                                  // the following line triggered an
+                                  // exception prior to r25670
+  tria.prepare_coarsening_and_refinement();
+  soltrans.prepare_for_coarsening_and_refinement(old_solution);
+  tria.execute_coarsening_and_refinement();
+  dof_handler.distribute_dofs (fe);
+  solution.reinit(dof_handler.n_dofs());
+  soltrans.interpolate(old_solution, solution);
+}
+
+
+int main()
+{
+  std::ofstream logfile("solution_transfer_02/output");
+  deallog.attach(logfile);
+  deallog.depth_console(0);
+  deallog.threshold_double(1.e-10);
+
+  deallog << "   1D solution transfer" << std::endl;
+  transfer<1>(logfile);
+
+  deallog << "   2D solution transfer" << std::endl;
+  transfer<2>(logfile);
+
+  deallog << "   3D solution transfer" << std::endl;
+  transfer<3>(logfile);
+}
+
+
+
diff --git a/tests/hp/solution_transfer_02/cmp/generic b/tests/hp/solution_transfer_02/cmp/generic
new file mode 100644 (file)
index 0000000..76a5bc7
--- /dev/null
@@ -0,0 +1,4 @@
+
+DEAL::   1D solution transfer
+DEAL::   2D solution transfer
+DEAL::   3D solution transfer

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.