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> > &)
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);
+ }
}
}
{
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);
--- /dev/null
+//---------------------------- 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);
+}
+
+
+
projects / dealii-svn.git / commitdiff