From: bangerth Date: Sat, 25 Jan 2014 04:11:46 +0000 (+0000) Subject: Clean up a good amount of code. In particular, do not try to access fe_indices of... X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=094a731e0e5774554fe45e287f4285c06e7d25f6;p=dealii-svn.git Clean up a good amount of code. In particular, do not try to access fe_indices of cells that are not active, since no finite element space is associated with them. git-svn-id: https://svn.dealii.org/trunk@32310 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/source/numerics/solution_transfer.cc b/deal.II/source/numerics/solution_transfer.cc index 54433972bb..d15d6768dd 100644 --- a/deal.II/source/numerics/solution_transfer.cc +++ b/deal.II/source/numerics/solution_transfer.cc @@ -234,9 +234,9 @@ void SolutionTransfer:: prepare_for_coarsening_and_refinement(const std::vector &all_in) { - Assert(prepared_for!=pure_refinement, ExcAlreadyPrepForRef()); - Assert(prepared_for!=coarsening_and_refinement, - ExcAlreadyPrepForCoarseAndRef()); + Assert (prepared_for!=pure_refinement, ExcAlreadyPrepForRef()); + Assert (prepared_for!=coarsening_and_refinement, + ExcAlreadyPrepForCoarseAndRef()); const unsigned int in_size=all_in.size(); Assert(in_size!=0, ExcNoInVectorsGiven()); @@ -244,7 +244,7 @@ prepare_for_coarsening_and_refinement(const std::vector &all_in) clear(); const unsigned int n_active_cells = dof_handler->get_tria().n_active_cells(); - n_dofs_old=dof_handler->n_dofs(); + n_dofs_old = dof_handler->n_dofs(); for (unsigned int i=0; i &all_in) } // first count the number - // of cells that'll be coarsened + // of cells that will be coarsened // and that'll stay or be refined unsigned int n_cells_to_coarsen=0; unsigned int n_cells_to_stay_or_refine=0; - typename DH::active_cell_iterator - act_cell = dof_handler->begin_active(), - endc = dof_handler->end(); - for (; act_cell!=endc; ++act_cell) + for (typename DH::active_cell_iterator act_cell = dof_handler->begin_active(); + act_cell!=dof_handler->end(); ++act_cell) { if (act_cell->coarsen_flag_set()) ++n_cells_to_coarsen; @@ -271,20 +269,16 @@ prepare_for_coarsening_and_refinement(const std::vector &all_in) ExcInternalError()); unsigned int n_coarsen_fathers=0; - typename DH::cell_iterator - cell=dof_handler->begin(); - for (; cell!=endc; ++cell) + for (typename DH::cell_iterator cell=dof_handler->begin(); + cell!=dof_handler->end(); ++cell) if (!cell->active() && cell->child(0)->coarsen_flag_set()) ++n_coarsen_fathers; - - if (n_cells_to_coarsen) - Assert(n_cells_to_coarsen>=2*n_coarsen_fathers, ExcInternalError()); + Assert(n_cells_to_coarsen>=2*n_coarsen_fathers, ExcInternalError()); // allocate the needed memory. initialize // the following arrays in an efficient // way, without copying much - std::vector > - (n_cells_to_stay_or_refine) + std::vector >(n_cells_to_stay_or_refine) .swap(indices_on_cell); std::vector > > @@ -293,19 +287,21 @@ prepare_for_coarsening_and_refinement(const std::vector &all_in) .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', unsigned int n_sr=0, n_cf=0; - cell = dof_handler->begin(); - for (; cell!=endc; ++cell) + for (typename DH::cell_iterator cell=dof_handler->begin(); + cell!=dof_handler->end(); ++cell) { + // CASE 1: active cell that remains as it is if (cell->active() && !cell->coarsen_flag_set()) { const unsigned int dofs_per_cell=cell->get_fe().dofs_per_cell; @@ -316,9 +312,11 @@ prepare_for_coarsening_and_refinement(const std::vector &all_in) // interpolating to the children cell->get_dof_indices(indices_on_cell[n_sr]); cell_map[std::make_pair(cell->level(), cell->index())] - = Pointerstruct(&indices_on_cell[n_sr],cell->active_fe_index()); + = Pointerstruct(&indices_on_cell[n_sr], cell->active_fe_index()); ++n_sr; } + + // CASE 2: cell is inactive but will become active else if (cell->has_children() && cell->child(0)->coarsen_flag_set()) { // note that if one child has the @@ -327,127 +325,51 @@ prepare_for_coarsening_and_refinement(const std::vector &all_in) // worked correctly for (unsigned int i=1; in_children(); ++i) Assert(cell->child(i)->coarsen_flag_set(), - ExcTriaPrepCoarseningNotCalledBefore()); - const unsigned int dofs_per_cell=cell->get_fe().dofs_per_cell; - - std::vector >(in_size, - Vector(dofs_per_cell)) - .swap(dof_values_on_cell[n_cf]); + ExcMessage("It looks like you didn't call " + "Triangulation::prepare_coarsening_and_refinement before " + "calling the current function. This can't work.")); + + // we will need to interpolate from the children of this cell + // to the current one. in the hp context, this also means + // we need to figure out which finite element space to interpolate + // to since that is not implied by the global FE as in the non-hp + // case. + bool different_fe_on_children = false; + for (unsigned int child=1; childn_children(); ++child) + if (cell->child(child)->active_fe_index() + != cell->child(0)->active_fe_index()) + { + different_fe_on_children = true; + break; + } - unsigned int fe_index = cell->active_fe_index(); + // take FE index from the child with most + // degrees of freedom locally 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 (different_fe_on_children == true) + for (unsigned int child=1; childn_children(); ++child) 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(); + const unsigned int target_fe_index = cell->child(most_general_child)->active_fe_index(); + + const unsigned int dofs_per_cell=cell->get_dof_handler().get_fe()[target_fe_index].dofs_per_cell; + std::vector >(in_size, + Vector(dofs_per_cell)) + .swap(dof_values_on_cell[n_cf]); + + + // store the data of each of the input vectors. get this data + // as interpolated onto a finite element space that encompasses + // that of all the children. note that cell->get_interpolated_dof_values + // already does all of the interpolations between spaces for (unsigned int j=0; jget_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); - - // 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); - - // now do the interpolation - // operation. FullMatrix only wants us - // to call vmult if the matrix size is - // actually non-zero, so check that - // case - const unsigned int dofs_per_cell = tmp2.size(); - tmp.reinit (dofs_per_cell, true); - if (dofs_per_cell > 0) - cell->child(most_general_child)->get_fe(). - get_restriction_matrix(child, cell->refinement_case()).vmult (tmp, tmp2); - for (unsigned int i=0; iget_interpolated_dof_values(all_in[j], + dof_values_on_cell[n_cf][j], + target_fe_index); cell_map[std::make_pair(cell->level(), cell->index())] - = Pointerstruct(&dof_values_on_cell[n_cf], fe_index); + = Pointerstruct(&dof_values_on_cell[n_cf], target_fe_index); ++n_cf; } }