SolutionTransfer<dim, VECTOR, DH>::
prepare_for_coarsening_and_refinement(const std::vector<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());
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<in_size; ++i)
{
}
// 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;
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<std::vector<types::global_dof_index> >
- (n_cells_to_stay_or_refine)
+ std::vector<std::vector<types::global_dof_index> >(n_cells_to_stay_or_refine)
.swap(indices_on_cell);
std::vector<std::vector<Vector<typename VECTOR::value_type> > >
.swap(dof_values_on_cell);
typename VECTOR::value_type zero_val = typename VECTOR::value_type();
+
Table<2,FullMatrix<double> > interpolation_hp;
std::vector<std::vector<bool> > restriction_is_additive;
+
internal::extract_interpolation_matrices (*dof_handler, interpolation_hp);
internal::restriction_additive (dof_handler->get_fe(), restriction_is_additive);
- Vector<typename VECTOR::value_type> 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;
// 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
// worked correctly
for (unsigned int i=1; i<cell->n_children(); ++i)
Assert(cell->child(i)->coarsen_flag_set(),
- ExcTriaPrepCoarseningNotCalledBefore());
- const unsigned int dofs_per_cell=cell->get_fe().dofs_per_cell;
-
- std::vector<Vector<typename VECTOR::value_type> >(in_size,
- Vector<typename VECTOR::value_type>(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; child<cell->n_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; child<cell->n_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; child<cell->n_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<Vector<typename VECTOR::value_type> >(in_size,
+ Vector<typename VECTOR::value_type>(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; j<in_size; ++j)
- {
- // store the data of each of
- // the input vectors
- 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; child<cell->n_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; i<dofs_per_cell; ++i)
- if (restriction_is_additive[fe_ind_general][i])
- dof_values_on_cell[n_cf][j](i) += tmp(i);
- else if (tmp(i) != zero_val)
- dof_values_on_cell[n_cf][j](i) = tmp(i);
- }
- }
- }
+ cell->get_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;
}
}
projects / dealii-svn.git / commitdiff