From 7b606f7c00e948ee089f384b2ae178cfbc4aafda Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Wed, 23 Sep 2015 13:28:13 -0500 Subject: [PATCH] Fix a place where hp::SolutionTransfer gets confused with FE_Nothing. This fixes a testcase reported by Claire Br. on the mailing list. --- doc/news/changes.h | 6 + source/dofs/dof_accessor_get.cc | 100 +++++++------ tests/hp/solution_transfer_13.cc | 210 +++++++++++++++++++++++++++ tests/hp/solution_transfer_13.output | 148 +++++++++++++++++++ 4 files changed, 418 insertions(+), 46 deletions(-) create mode 100644 tests/hp/solution_transfer_13.cc create mode 100644 tests/hp/solution_transfer_13.output diff --git a/doc/news/changes.h b/doc/news/changes.h index b96327e3f3..673990896e 100644 --- a/doc/news/changes.h +++ b/doc/news/changes.h @@ -356,6 +356,12 @@ inconvenience this causes.
    +
  1. Fixed: hp::SolutionTransfer could get confused when dealing with + FE_Nothing elements. This is now fixed. +
    + (Claire Br, Wolfgang Bangerth, 2015/09/23) +
    2. +
  2. Improved: The construction of the non-local graph for quick data exchange of TrilinosWrappers::SparseMatrix became very slow for a few thousand processors. This has been fixed. diff --git a/source/dofs/dof_accessor_get.cc b/source/dofs/dof_accessor_get.cc index ac4bcec477..9d248a8681 100644 --- a/source/dofs/dof_accessor_get.cc +++ b/source/dofs/dof_accessor_get.cc @@ -109,53 +109,61 @@ get_interpolated_dof_values (const InputVector &values, typename BaseClass::ExcVectorDoesNotMatch()); - Vector tmp1(dofs_per_cell); - Vector tmp2(dofs_per_cell); - - interpolated_values = 0; - - // later on we will have to push the values interpolated from the child - // to the mother cell into the output vector. unfortunately, there are - // two types of elements: ones where you add up the contributions from - // the different child cells, and ones where you overwrite. - // - // an example for the first is piecewise constant (and discontinuous) - // elements, where we build the value on the coarse cell by averaging - // the values from the cell (i.e. by adding up a fraction of the values - // of their values) - // - // an example for the latter are the usual continuous elements. the - // value on a vertex of a coarse cell must there be the same, - // irrespective of the adjacent cell we are presently on. so we always - // overwrite. in fact, we must, since we cannot know in advance how many - // neighbors there will be, so there is no way to compute the average - // with fixed factors - // - // so we have to find out to which type this element belongs. the - // difficulty is: the finite element may be a composed one, so we can - // only hope to do this for each shape function individually. in fact, - // there are even weird finite elements (for example the Raviart-Thomas - // element) which have shape functions that are additive (interior ones) - // and others that are overwriting (face degrees of freedom that need to - // be continuous across the face). - for (unsigned int child=0; childn_children(); ++child) + // see if the finite element we have on the current cell has any + // degrees of freedom to begin with; if not (e.g., when + // interpolating FE_Nothing), then simply skip all of the + // following since the output vector would be of size zero + // anyway (and in fact is of size zero, see the assertion above) + if (fe.dofs_per_cell > 0) { - // get the values from the present child, if necessary by - // interpolation itself either from its own children or - // by interpolating from the finite element on an active - // child to the finite element space requested here - this->child(child)->get_interpolated_dof_values (values, - tmp1, - fe_index); - // interpolate these to the mother cell - fe.get_restriction_matrix(child, this->refinement_case()).vmult (tmp2, tmp1); - - // and add up or set them in the output vector - for (unsigned int i=0; i tmp1(dofs_per_cell); + Vector tmp2(dofs_per_cell); + + interpolated_values = 0; + + // later on we will have to push the values interpolated from the child + // to the mother cell into the output vector. unfortunately, there are + // two types of elements: ones where you add up the contributions from + // the different child cells, and ones where you overwrite. + // + // an example for the first is piecewise constant (and discontinuous) + // elements, where we build the value on the coarse cell by averaging + // the values from the cell (i.e. by adding up a fraction of the values + // of their values) + // + // an example for the latter are the usual continuous elements. the + // value on a vertex of a coarse cell must there be the same, + // irrespective of the adjacent cell we are presently on. so we always + // overwrite. in fact, we must, since we cannot know in advance how many + // neighbors there will be, so there is no way to compute the average + // with fixed factors + // + // so we have to find out to which type this element belongs. the + // difficulty is: the finite element may be a composed one, so we can + // only hope to do this for each shape function individually. in fact, + // there are even weird finite elements (for example the Raviart-Thomas + // element) which have shape functions that are additive (interior ones) + // and others that are overwriting (face degrees of freedom that need to + // be continuous across the face). + for (unsigned int child=0; childn_children(); ++child) + { + // get the values from the present child, if necessary by + // interpolation itself either from its own children or + // by interpolating from the finite element on an active + // child to the finite element space requested here + this->child(child)->get_interpolated_dof_values (values, + tmp1, + fe_index); + // interpolate these to the mother cell + fe.get_restriction_matrix(child, this->refinement_case()).vmult (tmp2, tmp1); + + // and add up or set them in the output vector + for (unsigned int i=0; i + +#include + +#include +#include +#include +#include + +#include +#include +#include + +#include + +#include + +#include +#include +#include + +using namespace dealii; + +int main() +{ + std::ofstream logfile("output"); + deallog.attach(logfile); + deallog.depth_console(0); + deallog.threshold_double(1.e-10); + + Triangulation<2> triangulation(Triangulation<2>::none); + GridGenerator::hyper_cube (triangulation); + triangulation.refine_global(1); + + hp::FECollection<2> fe_collection; + fe_collection.push_back(FE_Q<2>(1)); + //fe_collection.push_back(FE_Q<2>(2)); + fe_collection.push_back(FE_Nothing<2>()); + + hp::DoFHandler<2> dof_handler(triangulation); + + // Assign FE + hp::DoFHandler<2>::active_cell_iterator cell = dof_handler.begin_active(); + hp::DoFHandler<2>::active_cell_iterator endc = dof_handler.end(); + + /* + * ----------- + * | 0 | 0 | + * ----------- + * | 1 | 1 | 0 - FEQ, 1 - FE_Nothing + * ----------- + */ + + cell->set_active_fe_index(1); + cell++; + cell->set_active_fe_index(1); + cell++; + cell->set_active_fe_index(0); + cell++; + cell->set_active_fe_index(0); + + dof_handler.distribute_dofs (fe_collection); + + // Init solution + Vector solution(dof_handler.n_dofs()); + solution = 1.0; + + + // Vector to visualize the FE of each cell + Vector FE_Type (triangulation.n_active_cells()); + unsigned int cnt_cells (0); + cell = dof_handler.begin_active(), + endc = dof_handler.end(); + for (; cell!=endc; ++cell) + { + unsigned int fe_index = cell->active_fe_index(); + FE_Type[cnt_cells] = fe_index; + ++ cnt_cells; + } + + // Save output + DataOut<2, hp::DoFHandler<2> > data_out; + data_out.attach_dof_handler (dof_handler); + data_out.add_data_vector (solution, "Solution"); + data_out.add_data_vector (FE_Type, "FE_Type"); + data_out.build_patches(); + data_out.write_gnuplot (deallog.get_file_stream()); + + + /* Set refine flags: + * ----------- + * | | R | + * ----------- + * | R | | + * ----------- + */ + + + cell = dof_handler.begin_active(); + cell->set_refine_flag(); + cell++; + cell++; + cell++; + cell->set_refine_flag(); + + triangulation.prepare_coarsening_and_refinement(); + + // Interpolate solution + SolutionTransfer<2, Vector, hp::DoFHandler<2> > solution_trans(dof_handler); + solution_trans.prepare_for_coarsening_and_refinement(solution); + + triangulation.execute_coarsening_and_refinement (); + + dof_handler.distribute_dofs(fe_collection); + + Vector new_solution(dof_handler.n_dofs()); + solution_trans.interpolate(solution, new_solution); + + FE_Type.reinit(triangulation.n_active_cells()); + cnt_cells=0; + cell = dof_handler.begin_active(), + endc = dof_handler.end(); + for (; cell!=endc; ++cell) + { + unsigned int fe_index = cell->active_fe_index(); + FE_Type[cnt_cells] = fe_index; + ++ cnt_cells; + } + + // Save new solution + DataOut<2, hp::DoFHandler<2> > data_out2; + data_out2.attach_dof_handler (dof_handler); + data_out2.add_data_vector (new_solution, "Solution"); + data_out2.add_data_vector(FE_Type, "FE_type" ); + data_out2.build_patches(); + data_out2.write_gnuplot(deallog.get_file_stream()); + + // Solution reinitialization + + dof_handler.distribute_dofs (fe_collection); + solution.reinit(dof_handler.n_dofs()); + solution = 1.0; + + /* Set coarsen flags: + * ----------- + * | | C | + * ----------- + * | C | | + * ----------- + */ + + endc = dof_handler.end(); + for (cell=dof_handler.begin_active(); cell!=endc; cell++) + if (cell->level() > 1) + cell->set_coarsen_flag(); + + triangulation.prepare_coarsening_and_refinement(); + + // Interpolate solution + SolutionTransfer<2, Vector, hp::DoFHandler<2> > solution_trans2(dof_handler); + solution_trans2.prepare_for_coarsening_and_refinement(solution); + + triangulation.execute_coarsening_and_refinement (); + + dof_handler.distribute_dofs (fe_collection); + + Vector new_solution2(dof_handler.n_dofs()); + solution_trans2.interpolate(solution, new_solution2); + + FE_Type.reinit(triangulation.n_active_cells()); + cnt_cells=0; + cell = dof_handler.begin_active(), + endc = dof_handler.end(); + for (; cell!=endc; ++cell) + { + unsigned int fe_index = cell->active_fe_index(); + FE_Type[cnt_cells] = fe_index; + ++ cnt_cells; + } + + // Save new solution + DataOut<2, hp::DoFHandler<2> > data_out3; + data_out3.attach_dof_handler (dof_handler); + data_out3.add_data_vector (new_solution2, "Solution"); + data_out3.add_data_vector(FE_Type, "FE_type" ); + data_out3.build_patches(); + data_out3.write_gnuplot(deallog.get_file_stream()); +} + diff --git a/tests/hp/solution_transfer_13.output b/tests/hp/solution_transfer_13.output new file mode 100644 index 0000000000..8994274ebc --- /dev/null +++ b/tests/hp/solution_transfer_13.output @@ -0,0 +1,148 @@ + +# This file was generated by the deal.II library. + + +# +# For a description of the GNUPLOT format see the GNUPLOT manual. +# +# +0.00000 0.00000 0.00000 1.00000 +0.500000 0.00000 0.00000 1.00000 + +0.00000 0.500000 0.00000 1.00000 +0.500000 0.500000 0.00000 1.00000 + + +0.500000 0.00000 0.00000 1.00000 +1.00000 0.00000 0.00000 1.00000 + +0.500000 0.500000 0.00000 1.00000 +1.00000 0.500000 0.00000 1.00000 + + +0.00000 0.500000 1.00000 0.00000 +0.500000 0.500000 1.00000 0.00000 + +0.00000 1.00000 1.00000 0.00000 +0.500000 1.00000 1.00000 0.00000 + + +0.500000 0.500000 1.00000 0.00000 +1.00000 0.500000 1.00000 0.00000 + +0.500000 1.00000 1.00000 0.00000 +1.00000 1.00000 1.00000 0.00000 + + +# This file was generated by the deal.II library. + + +# +# For a description of the GNUPLOT format see the GNUPLOT manual. +# +# +0.500000 0.00000 0.00000 1.00000 +1.00000 0.00000 0.00000 1.00000 + +0.500000 0.500000 0.00000 1.00000 +1.00000 0.500000 0.00000 1.00000 + + +0.00000 0.500000 1.00000 0.00000 +0.500000 0.500000 1.00000 0.00000 + +0.00000 1.00000 1.00000 0.00000 +0.500000 1.00000 1.00000 0.00000 + + +0.00000 0.00000 0.00000 1.00000 +0.250000 0.00000 0.00000 1.00000 + +0.00000 0.250000 0.00000 1.00000 +0.250000 0.250000 0.00000 1.00000 + + +0.250000 0.00000 0.00000 1.00000 +0.500000 0.00000 0.00000 1.00000 + +0.250000 0.250000 0.00000 1.00000 +0.500000 0.250000 0.00000 1.00000 + + +0.00000 0.250000 0.00000 1.00000 +0.250000 0.250000 0.00000 1.00000 + +0.00000 0.500000 0.00000 1.00000 +0.250000 0.500000 0.00000 1.00000 + + +0.250000 0.250000 0.00000 1.00000 +0.500000 0.250000 0.00000 1.00000 + +0.250000 0.500000 0.00000 1.00000 +0.500000 0.500000 0.00000 1.00000 + + +0.500000 0.500000 1.00000 0.00000 +0.750000 0.500000 1.00000 0.00000 + +0.500000 0.750000 1.00000 0.00000 +0.750000 0.750000 1.00000 0.00000 + + +0.750000 0.500000 1.00000 0.00000 +1.00000 0.500000 1.00000 0.00000 + +0.750000 0.750000 1.00000 0.00000 +1.00000 0.750000 1.00000 0.00000 + + +0.500000 0.750000 1.00000 0.00000 +0.750000 0.750000 1.00000 0.00000 + +0.500000 1.00000 1.00000 0.00000 +0.750000 1.00000 1.00000 0.00000 + + +0.750000 0.750000 1.00000 0.00000 +1.00000 0.750000 1.00000 0.00000 + +0.750000 1.00000 1.00000 0.00000 +1.00000 1.00000 1.00000 0.00000 + + +# This file was generated by the deal.II library. + + +# +# For a description of the GNUPLOT format see the GNUPLOT manual. +# +# +0.00000 0.00000 0.00000 1.00000 +0.500000 0.00000 0.00000 1.00000 + +0.00000 0.500000 0.00000 1.00000 +0.500000 0.500000 0.00000 1.00000 + + +0.500000 0.00000 0.00000 1.00000 +1.00000 0.00000 0.00000 1.00000 + +0.500000 0.500000 0.00000 1.00000 +1.00000 0.500000 0.00000 1.00000 + + +0.00000 0.500000 1.00000 0.00000 +0.500000 0.500000 1.00000 0.00000 + +0.00000 1.00000 1.00000 0.00000 +0.500000 1.00000 1.00000 0.00000 + + +0.500000 0.500000 1.00000 0.00000 +1.00000 0.500000 1.00000 0.00000 + +0.500000 1.00000 1.00000 0.00000 +1.00000 1.00000 1.00000 0.00000 + + -- 2.39.5