From 729e38bc471ea98ba0ef309bd81e711bdfafede8 Mon Sep 17 00:00:00 2001 From: kronbichler Date: Sat, 10 Mar 2012 11:22:23 +0000 Subject: [PATCH] Solution transfer for hp elements. git-svn-id: https://svn.dealii.org/trunk@25237 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/doc/news/changes.h | 6 + .../deal.II/numerics/solution_transfer.h | 18 +- deal.II/source/numerics/solution_transfer.cc | 251 +++++++++++++++--- tests/hp/solution_transfer.cc | 246 +++++++++++++++++ tests/hp/solution_transfer/cmp/generic | 10 + 5 files changed, 498 insertions(+), 33 deletions(-) create mode 100644 tests/hp/solution_transfer.cc create mode 100644 tests/hp/solution_transfer/cmp/generic diff --git a/deal.II/doc/news/changes.h b/deal.II/doc/news/changes.h index 03e8d52645..c16b1b6e5f 100644 --- a/deal.II/doc/news/changes.h +++ b/deal.II/doc/news/changes.h @@ -152,6 +152,12 @@ enabled due to a missing include file in file

Specific improvements

    +
  1. Extended: +SolutionTransfer is now also able to transfer solutions between hp::DoFHandler where +the finite element index changes during refinement. +
    +(Katharina Kormann, Martin Kronbichler, 2012/03/10) +
  2. Changed: A new method to determine an initial guess for the Newton method was coded in MappingQ::transform_real_to_unit_cell. diff --git a/deal.II/include/deal.II/numerics/solution_transfer.h b/deal.II/include/deal.II/numerics/solution_transfer.h index a214daaeb6..1e9a1c31f6 100644 --- a/deal.II/include/deal.II/numerics/solution_transfer.h +++ b/deal.II/include/deal.II/numerics/solution_transfer.h @@ -273,7 +273,7 @@ class SolutionTransfer /** * This function * interpolates the discrete functions - * that are stored in @p all_out onto + * that are stored in @p all_in onto * the refined and/or coarsenend grid. * It assumes the vectors in @p all_in * denote the same vectors @@ -423,17 +423,29 @@ class SolutionTransfer * vector dof_values (if the * children of this cell will be deleted) * is needed, hence one @p user_pointer should - * be sufficient, but to allow some errorchecks + * be sufficient, but to allow some error checks * and to preserve the user from making * user errors the @p user_pointer will be * 'multiplied' by this structure. */ struct Pointerstruct { - Pointerstruct(); + Pointerstruct() : indices_ptr(0), dof_values_ptr(0), active_fe_index(0) {}; + Pointerstruct(std::vector *indices_ptr_in, + const unsigned int active_fe_index_in = 0) + : + indices_ptr(indices_ptr_in), + dof_values_ptr (0), + active_fe_index(active_fe_index_in) {}; + Pointerstruct(std::vector > *dof_values_ptr_in, + const unsigned int active_fe_index_in = 0) : + indices_ptr (0), + dof_values_ptr(dof_values_ptr_in), + active_fe_index(active_fe_index_in) {}; std::size_t memory_consumption () const; std::vector *indices_ptr; std::vector > *dof_values_ptr; + unsigned int active_fe_index; }; /** diff --git a/deal.II/source/numerics/solution_transfer.cc b/deal.II/source/numerics/solution_transfer.cc index eeca9f945e..a88ae40eda 100644 --- a/deal.II/source/numerics/solution_transfer.cc +++ b/deal.II/source/numerics/solution_transfer.cc @@ -40,6 +40,7 @@ SolutionTransfer::SolutionTransfer(const DH &dof) {} + template SolutionTransfer::~SolutionTransfer() { @@ -47,13 +48,6 @@ SolutionTransfer::~SolutionTransfer() } -template -SolutionTransfer::Pointerstruct::Pointerstruct(): - indices_ptr(0), - dof_values_ptr(0) -{} - - template void SolutionTransfer::clear () @@ -66,6 +60,7 @@ void SolutionTransfer::clear () } + template void SolutionTransfer::prepare_for_pure_refinement() { @@ -101,6 +96,7 @@ void SolutionTransfer::prepare_for_pure_refinement() } + template void SolutionTransfer::refine_interpolate(const VECTOR &in, @@ -155,6 +151,51 @@ SolutionTransfer::refine_interpolate(const VECTOR &in, } + +namespace internal +{ + template + void extract_interpolation_matrices (const DH&, + Table<2,FullMatrix > &) + {} + + template + void extract_interpolation_matrices (const dealii::hp::DoFHandler &dof, + Table<2,FullMatrix > &matrices) + { + const dealii::hp::FECollection &fe = dof.get_fe(); + matrices.reinit (fe.size(), fe.size()); + for (unsigned int i=0; i + void restriction_additive (const FiniteElement &, + std::vector > &) + {} + + template + void restriction_additive (const dealii::hp::FECollection &fe, + std::vector > &restriction_is_additive) + { + restriction_is_additive.resize (fe.size()); + for (unsigned int f=0; f void SolutionTransfer:: @@ -218,6 +259,13 @@ prepare_for_coarsening_and_refinement(const std::vector &all_in) std::vector > (in_size)) .swap(dof_values_on_cell); + typename VECTOR::value_type zero_val = typename VECTOR::value_type(); + Table<2,FullMatrix > interpolation_hp; + std::vector > restriction_is_additive; + internal::extract_interpolation_matrices (*dof_handler, interpolation_hp); + internal::restriction_additive (dof_handler->get_fe(), restriction_is_additive); + Vector tmp, tmp2; + // we need counters for // the 'to_stay_or_refine' cells 'n_sr' and // the 'coarsen_fathers' cells 'n_cf', @@ -234,7 +282,8 @@ prepare_for_coarsening_and_refinement(const std::vector &all_in) // dof indices and later // interpolating to the children cell->get_dof_indices(indices_on_cell[n_sr]); - cell_map[std::make_pair(cell->level(), cell->index())].indices_ptr=&indices_on_cell[n_sr]; + cell_map[std::make_pair(cell->level(), cell->index())] + = Pointerstruct(&indices_on_cell[n_sr],cell->active_fe_index()); ++n_sr; } else if (cell->has_children() && cell->child(0)->coarsen_flag_set()) @@ -249,17 +298,76 @@ prepare_for_coarsening_and_refinement(const std::vector &all_in) const unsigned int dofs_per_cell=cell->get_fe().dofs_per_cell; std::vector >(in_size, - Vector(dofs_per_cell)) + Vector(dofs_per_cell)) .swap(dof_values_on_cell[n_cf]); + unsigned int fe_index = cell->active_fe_index(); + unsigned int most_general_child = 0; + bool different_elements = false; + for (unsigned int child=0; childn_children(); ++child) + { + if (cell->child(child)->active_fe_index() != fe_index) + different_elements = true; + // take FE index from the child with most + // degrees of freedom locally + if (cell->child(child)->get_fe().dofs_per_cell > + cell->child(most_general_child)->get_fe().dofs_per_cell) + most_general_child = child; + } + if (different_elements == true) + fe_index = cell->child(most_general_child)->active_fe_index(); + for (unsigned int j=0; jget_interpolated_dof_values(all_in[j], - dof_values_on_cell[n_cf][j]); + if (different_elements == false) + cell->get_interpolated_dof_values(all_in[j], + dof_values_on_cell[n_cf][j]); + else + { + // if we have different elements, first + // interpolate the children's contribution to + // the most general FE on the child + // level. Then we manually write the + // interpolation operation to the coarser + // level + dof_values_on_cell[n_cf][j].reinit (cell->child(most_general_child)->get_fe().dofs_per_cell); + const unsigned int fe_ind_general = + cell->child(most_general_child)->active_fe_index(); + for (unsigned int child=0; childn_children(); ++child) + { + tmp.reinit (cell->child(child)->get_fe().dofs_per_cell, + true); + cell->child(child)->get_dof_values (all_in[j], + tmp); + const unsigned int child_ind = + cell->child(child)->active_fe_index(); + if (child_ind != fe_ind_general) + { + tmp2.reinit (cell->child(most_general_child)->get_fe().dofs_per_cell, + true); + interpolation_hp (fe_ind_general, child_ind).vmult (tmp2, tmp); + } + else + tmp2.swap (tmp); + + // now do the interpolation operation + const unsigned int dofs_per_cell = tmp2.size(); + tmp.reinit (dofs_per_cell, true); + cell->child(most_general_child)->get_fe(). + get_restriction_matrix(child, cell->refinement_case()).vmult (tmp, tmp2); + for (unsigned int i=0; ilevel(), cell->index())].dof_values_ptr=&dof_values_on_cell[n_cf]; + cell_map[std::make_pair(cell->level(), cell->index())] + = Pointerstruct(&dof_values_on_cell[n_cf], fe_index); ++n_cf; } } @@ -270,6 +378,7 @@ prepare_for_coarsening_and_refinement(const std::vector &all_in) } + template void SolutionTransfer::prepare_for_coarsening_and_refinement(const VECTOR &in) @@ -279,6 +388,7 @@ SolutionTransfer::prepare_for_coarsening_and_refinement(const V } + template void SolutionTransfer:: interpolate (const std::vector &all_in, @@ -306,6 +416,10 @@ interpolate (const std::vector &all_in, pointerstruct, cell_map_end=cell_map.end(); + Table<2,FullMatrix > interpolation_hp; + internal::extract_interpolation_matrices (*dof_handler, interpolation_hp); + Vector tmp, tmp2; + typename DH::cell_iterator cell = dof_handler->begin(), endc = dof_handler->end(); for (; cell!=endc; ++cell) @@ -328,28 +442,91 @@ interpolate (const std::vector &all_in, { Assert (valuesptr == 0, ExcInternalError()); - // make sure that the size of the - // stored indices is the same as - // dofs_per_cell. this is kind of - // a test if we use the same fe - // in the hp case. to really do - // that test we would have to - // store the fe_index of all - // cells - Assert(dofs_per_cell==(*indexptr).size(), - ExcNumberOfDoFsPerCellHasChanged()); + // get the values of // each of the input // data vectors on this // cell and prolong it // to its children + unsigned int in_size = indexptr->size(); local_values.reinit(dofs_per_cell, true); - for (unsigned int j=0; jset_dof_values_by_interpolation(local_values, - all_out[j]); + // check the FE index of the new element. if + // we have children and not all of the + // children have the same FE index, need to + // manually implement + // set_dof_values_by_interpolation + unsigned int new_fe_index = cell->active_fe_index(); + bool different_elements = false; + if (cell->has_children()) + { + new_fe_index = cell->child(0)->active_fe_index(); + for (unsigned int child=1; childn_children(); ++child) + if (cell->child(child)->active_fe_index() != new_fe_index) + { + different_elements = true; + break; + } + } + if (different_elements == true || + new_fe_index != pointerstruct->second.active_fe_index) + { + const unsigned int old_index = + pointerstruct->second.active_fe_index; + tmp.reinit (in_size, true); + for (unsigned int i=0; ihas_children() ? + cell->child(0)->get_fe().dofs_per_cell + : cell->get_fe().dofs_per_cell, true); + AssertDimension (local_values.size(), + interpolation_hp(new_fe_index,old_index).m()); + // simple case where all children have the + // same FE index: just interpolate to their FE + // first and then use the standard routines + interpolation_hp(new_fe_index,old_index).vmult (local_values, tmp); + } + + if (cell->has_children() == false) + cell->set_dof_values (local_values, all_out[j]); + else + for (unsigned int child=0; childn_children(); ++child) + { + if (different_elements == true) + { + const unsigned int c_index = + cell->child(child)->active_fe_index(); + if (c_index != old_index) + { + AssertDimension (tmp.size(), + interpolation_hp(c_index,old_index).n()); + local_values.reinit(cell->child(child)->get_fe().dofs_per_cell, true); + AssertDimension (local_values.size(), + interpolation_hp(c_index,old_index).m()); + interpolation_hp(c_index,old_index).vmult (local_values, tmp); + } + else + local_values = tmp; + } + tmp2.reinit (cell->child(child)->get_fe().dofs_per_cell, true); + cell->child(child)->get_fe().get_prolongation_matrix(child, cell->refinement_case()) + .vmult (tmp2, local_values); + cell->child(child)->set_dof_values(tmp2,all_out[j]); + } + } + else + { + AssertDimension (dofs_per_cell, indexptr->size()); + for (unsigned int i=0; iset_dof_values_by_interpolation(local_values, + all_out[j]); + } } } // children of cell were @@ -379,11 +556,25 @@ interpolate (const std::vector &all_in, // test we would have to // store the fe_index of all // cells - Assert(dofs_per_cell==(*valuesptr)[j].size(), - ExcNumberOfDoFsPerCellHasChanged()); + const Vector* data = 0; + const unsigned int active_fe_index = cell->active_fe_index(); + if (active_fe_index != pointerstruct->second.active_fe_index) + { + const unsigned int old_index = pointerstruct->second.active_fe_index; + tmp.reinit (dofs_per_cell, true); + AssertDimension ((*valuesptr)[j].size(), + interpolation_hp(active_fe_index,old_index).n()); + AssertDimension (tmp.size(), + interpolation_hp(active_fe_index,old_index).m()); + interpolation_hp(active_fe_index,old_index).vmult (tmp, (*valuesptr)[j]); + data = &tmp; + } + else + data = &(*valuesptr)[j]; + for (unsigned int i=0; i +#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 +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 function; + Triangulation tria; + GridGenerator::hyper_cube(tria); + tria.refine_global(5-dim); + const unsigned int max_degree = 6-dim; + hp::FECollection fe_q; + hp::FECollection fe_dgq; + for (unsigned int deg=1; deg<=max_degree; ++deg) + { + fe_q.push_back (FE_Q(deg)); + fe_dgq.push_back (FE_DGQ(deg)); + } + hp::DoFHandler q_dof_handler(tria); + hp::DoFHandler dgq_dof_handler(tria); + Vector q_solution; + Vector dgq_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(), + celldg=dgq_dof_handler.begin_active(), + endc=q_dof_handler.end(); + for (; cell!=endc; ++cell, ++celldg, ++counter) + { + if (counter < 15) + cell->set_active_fe_index(1); + else + cell->set_active_fe_index(rand()%max_degree); + if (counter < 15) + celldg->set_active_fe_index(1); + else + celldg->set_active_fe_index(rand()%max_degree); + } + } + + q_dof_handler.distribute_dofs (fe_q); + q_solution.reinit(q_dof_handler.n_dofs()); + + dgq_dof_handler.distribute_dofs (fe_dgq); + dgq_solution.reinit(dgq_dof_handler.n_dofs()); + + ConstraintMatrix cm; + cm.close(); + VectorTools::interpolate (mapping, q_dof_handler, function, q_solution); + VectorTools::interpolate (mapping, dgq_dof_handler, function, dgq_solution); + + SolutionTransfer,hp::DoFHandler > q_soltrans(q_dof_handler); + SolutionTransfer,hp::DoFHandler > dgq_soltrans(dgq_dof_handler); + + + // test b): do some coarsening and + // refinement + q_soltrans.clear(); + dgq_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 (rand() % 3 == 0) + cell->set_refine_flag(); + else if (rand() % 3 == 3) + cell->set_coarsen_flag(); + } + + Vector q_old_solution=q_solution, + dgq_old_solution=dgq_solution; + tria.prepare_coarsening_and_refinement(); + q_soltrans.prepare_for_coarsening_and_refinement(q_old_solution); + dgq_soltrans.prepare_for_coarsening_and_refinement(dgq_old_solution); + tria.execute_coarsening_and_refinement(); + + counter = 0; + { + typename hp::DoFHandler::active_cell_iterator + cell = q_dof_handler.begin_active(), + celldg = dgq_dof_handler.begin_active(), + endc = q_dof_handler.end(); + for (; cell!=endc; ++cell, ++celldg, ++counter) + { + if (counter > 20 && counter < 90) + cell->set_active_fe_index(0); + else + cell->set_active_fe_index(rand()%max_degree); + if (counter > 20 && counter < 90) + celldg->set_active_fe_index(0); + else + celldg->set_active_fe_index(rand()%max_degree); + } + } + + q_dof_handler.distribute_dofs (fe_q); + dgq_dof_handler.distribute_dofs (fe_dgq); + q_solution.reinit(q_dof_handler.n_dofs()); + dgq_solution.reinit(dgq_dof_handler.n_dofs()); + q_soltrans.interpolate(q_old_solution, q_solution); + dgq_soltrans.interpolate(dgq_old_solution, dgq_solution); + + // check correctness by comparing the values + // on points of QGauss of order 2. + MyFunction func; + { + 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 quad (QGauss (2)); + hp::FEValues hp_fe_val (fe_dgq, quad, update_values | + update_quadrature_points); + std::vector vals (quad[0].size()); + typename hp::DoFHandler::active_cell_iterator + celldg = dgq_dof_handler.begin_active(), + endc = dgq_dof_handler.end(); + for (; celldg!=endc; ++celldg) + { + hp_fe_val.reinit (celldg, 0); + const FEValues &fe_val = hp_fe_val.get_present_fe_values(); + fe_val.get_function_values (dgq_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/cmp/generic b/tests/hp/solution_transfer/cmp/generic new file mode 100644 index 0000000000..7af228e047 --- /dev/null +++ b/tests/hp/solution_transfer/cmp/generic @@ -0,0 +1,10 @@ + +DEAL:: 1D solution transfer +DEAL::Error in interpolating hp FE_Q: 0 +DEAL::Error in interpolating hp FE_DGQ: 0 +DEAL:: 2D solution transfer +DEAL::Error in interpolating hp FE_Q: 0 +DEAL::Error in interpolating hp FE_DGQ: 0 +DEAL:: 3D solution transfer +DEAL::Error in interpolating hp FE_Q: 0 +DEAL::Error in interpolating hp FE_DGQ: 0 -- 2.39.5