From bba3007ca8f54d9bb88f6ff0813fcc188a178994 Mon Sep 17 00:00:00 2001 From: bangerth Date: Fri, 29 Jun 2012 08:25:48 +0000 Subject: [PATCH] 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. git-svn-id: https://svn.dealii.org/trunk@25660 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/doc/news/changes.h | 6 ++ deal.II/source/numerics/solution_transfer.cc | 86 ++++++++++++++- tests/hp/solution_transfer_02.cc | 108 +++++++++++++++++++ tests/hp/solution_transfer_02/cmp/generic | 4 + 4 files changed, 202 insertions(+), 2 deletions(-) create mode 100644 tests/hp/solution_transfer_02.cc create mode 100644 tests/hp/solution_transfer_02/cmp/generic diff --git a/deal.II/doc/news/changes.h b/deal.II/doc/news/changes.h index b309261733..49d651cd18 100644 --- a/deal.II/doc/news/changes.h +++ b/deal.II/doc/news/changes.h @@ -47,6 +47,12 @@ used for boundary indicators.
    +
  1. 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. +
    +(Wolfgang Bangerth, 2012/06/29) +
  2. 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 diff --git a/deal.II/source/numerics/solution_transfer.cc b/deal.II/source/numerics/solution_transfer.cc index cb6f005551..f779990306 100644 --- a/deal.II/source/numerics/solution_transfer.cc +++ b/deal.II/source/numerics/solution_transfer.cc @@ -154,6 +154,19 @@ SolutionTransfer::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 void extract_interpolation_matrices (const DH&, Table<2,FullMatrix > &) @@ -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::ExcInterpolationNotImplemented &) + { + matrices(i,j).reinit (0,0); + } } } @@ -347,7 +386,50 @@ prepare_for_coarsening_and_refinement(const std::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 index 0000000000..33de200516 --- /dev/null +++ b/tests/hp/solution_transfer_02.cc @@ -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 +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + + +template +void transfer(std::ostream &out) +{ + Triangulation tria; + GridGenerator::hyper_cube(tria); + tria.refine_global(1); + + hp::FECollection fe; + fe.push_back (FE_Q (1)); + fe.push_back (FE_Nothing()); + + // create a DoFHandler on which we + // have both cells with FE_Q as + // well as FE_Nothing + hp::DoFHandler dof_handler(tria); + dof_handler.begin(0)->child(0)->set_active_fe_index(1); + + Vector solution; + ConstraintMatrix cm; + cm.close(); + + dof_handler.distribute_dofs (fe); + solution.reinit(dof_handler.n_dofs()); + + for (unsigned int i=0; i,hp::DoFHandler > soltrans(dof_handler); + + typename Triangulation::active_cell_iterator cell=tria.begin_active(), + endc=tria.end(); + ++cell; + ++cell; + for (; cell!=endc; ++cell) + cell->set_refine_flag(); + + Vector 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 index 0000000000..76a5bc702c --- /dev/null +++ b/tests/hp/solution_transfer_02/cmp/generic @@ -0,0 +1,4 @@ + +DEAL:: 1D solution transfer +DEAL:: 2D solution transfer +DEAL:: 3D solution transfer -- 2.39.5