From 8f4f6198a0c421ec6559d8dbf255b6d6fb722e81 Mon Sep 17 00:00:00 2001 From: Denis Davydov Date: Mon, 13 Apr 2015 11:19:11 +0200 Subject: [PATCH] added a SolutionTransfer test for FE_Q_Hierarchical. Also cleaned up unnecessary Assert in the class. --- source/fe/fe_q_hierarchical.cc | 4 - tests/hp/solution_transfer_12.cc | 238 +++++++++++++++++++++++++++ tests/hp/solution_transfer_12.output | 7 + 3 files changed, 245 insertions(+), 4 deletions(-) create mode 100644 tests/hp/solution_transfer_12.cc create mode 100644 tests/hp/solution_transfer_12.output diff --git a/source/fe/fe_q_hierarchical.cc b/source/fe/fe_q_hierarchical.cc index d4a4017ee2..90e6d1ed71 100644 --- a/source/fe/fe_q_hierarchical.cc +++ b/source/fe/fe_q_hierarchical.cc @@ -164,10 +164,6 @@ FE_Q_Hierarchical::get_interpolation_matrix(const FiniteElement< dim> &sour ExcDimensionMismatch (matrix.m(), source_fe->dofs_per_cell)); - // make sure that we interpolate to a richer element. - Assert (this->dofs_per_cell >= source_fe->dofs_per_cell, - ExcMessage("Interpolation to higher polynomial degree is advised.")); - // Recall that DoFs are renumbered in the following order: // vertices, lines, quads, hexes. // As we deal with hierarchical FE, interpolation matrix is rather easy: diff --git a/tests/hp/solution_transfer_12.cc b/tests/hp/solution_transfer_12.cc new file mode 100644 index 0000000000..c2e29e3fda --- /dev/null +++ b/tests/hp/solution_transfer_12.cc @@ -0,0 +1,238 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 1998 - 2014 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. +// +// --------------------------------------------------------------------- + + + +#include "../tests.h" +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +// a linear function that should be transferred exactly with Q1 and Q2 +// elements +// This is a copy-paste of solution_transfer.cc test, +// with the two exceptions: +// - FE_Q_Hierarchical +// - manual interpolation as there is no general interpolation implemented +// for this FE currently, thus VectorTools::interpolate can't be used. +template +class MyFunction : public Function +{ +public: + MyFunction () : Function() {}; + + virtual double value (const Point &p, + const unsigned int) const + { + double f=0.25 + 2 * p[0]; + if (dim>1) + f+=0.772 * p[1]; + if (dim>2) + f-=3.112 * p[2]; + return f; + }; +}; + + +template +void transfer(std::ostream &out) +{ + MyFunction func; + Triangulation tria; + GridGenerator::hyper_cube(tria); + tria.refine_global(5-dim); + const unsigned int max_degree = 6-dim; + hp::FECollection fe_q; + for (unsigned int deg=1; deg<=max_degree; ++deg) + { + fe_q.push_back (FE_Q_Hierarchical(deg)); + } + hp::DoFHandler q_dof_handler(tria); + Vector q_solution; + MappingQ1 mapping; + + // refine a few cells + typename Triangulation::active_cell_iterator cell=tria.begin_active(), + endc=tria.end(); + ++cell; + ++cell; + for (; cell!=endc; ++cell) + cell->set_refine_flag(); + tria.prepare_coarsening_and_refinement(); + tria.execute_coarsening_and_refinement(); + + // randomly assign FE orders + unsigned int counter = 0; + { + typename hp::DoFHandler::active_cell_iterator + cell=q_dof_handler.begin_active(), + endc=q_dof_handler.end(); + for (; cell!=endc; ++cell, ++counter) + { + if (counter < 15) + cell->set_active_fe_index(1); + else + cell->set_active_fe_index(Testing::rand()%max_degree); + } + } + + q_dof_handler.distribute_dofs (fe_q); + q_solution.reinit(q_dof_handler.n_dofs()); + + ConstraintMatrix cm; + cm.close(); + + // interpolation is not straight forward for FE_Q_Hierarchical. + // Since we know that we do linear function, just get those DoFs + // which correspond to vertices and assign a value of shape function there + // while keeping other modes zero. + { + q_solution = 0.; + + hp::QCollection quad; + quad.push_back (QTrapez ()); + hp::FEValues hp_fe_val (fe_q, quad, update_values | + update_quadrature_points); + std::vector local_dof_indices; + typename hp::DoFHandler::active_cell_iterator + cell = q_dof_handler.begin_active(), + endc = q_dof_handler.end(); + for (; cell!=endc; ++cell) + { + hp_fe_val.reinit (cell, 0); + const FEValues &fe_val = hp_fe_val.get_present_fe_values(); + local_dof_indices.resize(fe_val.dofs_per_cell); + + cell->get_dof_indices(local_dof_indices); + + Assert (local_dof_indices.size() >= fe_val.n_quadrature_points, + ExcInternalError()); + + for (unsigned int q=0; q,hp::DoFHandler > q_soltrans(q_dof_handler); + + + // test b): do some coarsening and + // refinement + q_soltrans.clear(); + + counter = 0; + cell=tria.begin_active(); + endc=tria.end(); + for (; cell!=endc; ++cell, ++counter) + { + if (counter > 120) + cell->set_coarsen_flag(); + else if (Testing::rand() % 3 == 0) + cell->set_refine_flag(); + else if (Testing::rand() % 3 == 3) + cell->set_coarsen_flag(); + } + + Vector q_old_solution=q_solution; + tria.prepare_coarsening_and_refinement(); + q_soltrans.prepare_for_coarsening_and_refinement(q_old_solution); + tria.execute_coarsening_and_refinement(); + + counter = 0; + { + typename hp::DoFHandler::active_cell_iterator + cell = q_dof_handler.begin_active(), + endc = q_dof_handler.end(); + for (; cell!=endc; ++cell, ++counter) + { + if (counter > 20 && counter < 90) + cell->set_active_fe_index(0); + else + cell->set_active_fe_index(Testing::rand()%max_degree); + } + } + + q_dof_handler.distribute_dofs (fe_q); + q_solution.reinit(q_dof_handler.n_dofs()); + q_soltrans.interpolate(q_old_solution, q_solution); + + // check correctness by comparing the values + // on points of QGauss of order 2. + { + double error = 0; + const hp::QCollection quad (QGauss (2)); + hp::FEValues hp_fe_val (fe_q, quad, update_values | + update_quadrature_points); + std::vector vals (quad[0].size()); + typename hp::DoFHandler::active_cell_iterator + cell = q_dof_handler.begin_active(), + endc = q_dof_handler.end(); + for (; cell!=endc; ++cell) + { + hp_fe_val.reinit (cell, 0); + const FEValues &fe_val = hp_fe_val.get_present_fe_values(); + fe_val.get_function_values (q_solution, vals); + for (unsigned int q=0; q(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_12.output b/tests/hp/solution_transfer_12.output new file mode 100644 index 0000000000..adc30e4418 --- /dev/null +++ b/tests/hp/solution_transfer_12.output @@ -0,0 +1,7 @@ + +DEAL:: 1D solution transfer +DEAL::Error in interpolating hp FE_Q_Hierarchical: 0 +DEAL:: 2D solution transfer +DEAL::Error in interpolating hp FE_Q_Hierarchical: 0 +DEAL:: 3D solution transfer +DEAL::Error in interpolating hp FE_Q_Hierarchical: 0 -- 2.39.5