get_interpolation_matrix (const FiniteElement<dim,spacedim> &x_source_fe,
FullMatrix<double> &interpolation_matrix) const
{
- Assert (interpolation_matrix.m() == this->dofs_per_cell,
- ExcDimensionMismatch (interpolation_matrix.m(),
- this->dofs_per_cell));
- Assert (interpolation_matrix.n() == x_source_fe.dofs_per_cell,
- ExcDimensionMismatch (interpolation_matrix.m(),
- x_source_fe.dofs_per_cell));
-
// go through the list of elements we can interpolate from
if (const FE_Q_Base<POLY,dim,spacedim> *source_fe
= dynamic_cast<const FE_Q_Base<POLY,dim,spacedim>*>(&x_source_fe))
{
// ok, source is a Q element, so we will be able to do the work
+ Assert (interpolation_matrix.m() == this->dofs_per_cell,
+ ExcDimensionMismatch (interpolation_matrix.m(),
+ this->dofs_per_cell));
+ Assert (interpolation_matrix.n() == x_source_fe.dofs_per_cell,
+ ExcDimensionMismatch (interpolation_matrix.m(),
+ x_source_fe.dofs_per_cell));
// only evaluate Q dofs
const unsigned int q_dofs_per_cell = Utilities::fixed_power<dim>(this->degree+1);
// element represents a function that is constant zero and has no
// degrees of freedom, so the interpolation is simply a multiplication
// with a n_dofs x 0 matrix. there is nothing to do here
+
+ // we would like to verify that the number of rows and columns of
+ // the matrix equals this->dofs_per_cell and zero. unfortunately,
+ // whenever we do FullMatrix::reinit(m,0), it sets both rows and
+ // columns to zero, instead of m and zero. thus, only test the
+ // number of columns
+ Assert (interpolation_matrix.n() == x_source_fe.dofs_per_cell,
+ ExcDimensionMismatch (interpolation_matrix.m(),
+ x_source_fe.dofs_per_cell));
+
}
else
AssertThrow (false,
// back to size zero. so if we get here and that is
// the wrong size, then this may be because the elements
// haven't implemented the correct function yet
+ //
+ // there is one wrinkle. we would like to only error out if
+ // the size of the matrix is 0 times 0 but at least one
+ // of the elements has more than one dof per cell. the
+ // problem is that if you reinit a matrix to 4x0, it automatically
+ // sets its size to 0x0. so we can only execute the following
+ // test if *both* elements have dofs_per_cell>0, not if *at
+ // least one* have.
Assert (! ((interpolation_hp(new_fe_index,old_index).m() == 0)
&&
(interpolation_hp(new_fe_index,old_index).n() == 0)
&&
((dof_handler->get_fe()[new_fe_index].dofs_per_cell > 0)
- ||
+ &&
(dof_handler->get_fe()[old_index].dofs_per_cell > 0))),
ExcMessage ("The interpolation between two different "
"elements you are trying to use here has "
"not been implemented for this pair of "
- "elememts!"));
+ "elements!"));
- 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 (tmp.size() > 0)
+ interpolation_hp(new_fe_index,old_index).vmult (local_values, tmp);
+ else
+ local_values = 0;
}
if (cell->has_children() == false)
(interpolation_hp(c_index,old_index).n() == 0)
&&
((dof_handler->get_fe()[c_index].dofs_per_cell > 0)
- ||
+ &&
(dof_handler->get_fe()[old_index].dofs_per_cell > 0))),
ExcMessage ("The interpolation between two different "
"elements you are trying to use here has "
"not been implemented for this pair of "
- "elememts!"));
-
- AssertDimension (local_values.size(),
- interpolation_hp(c_index,old_index).m());
- interpolation_hp(c_index,old_index).vmult (local_values, tmp);
+ "elements!"));
+
+ if (tmp.size() > 0)
+ interpolation_hp(c_index,old_index).vmult (local_values, tmp);
+ else
+ local_values = 0;
}
else
local_values = tmp;
triangulation.refine_global (1);
- hp::DoFHandler<2> dof_handler(triangulation);
hp::FECollection<2> fe_collection;
-
-
fe_collection.push_back(FE_Q<2>(1));
fe_collection.push_back(FE_Nothing<2>());
+ hp::DoFHandler<2> dof_handler(triangulation);
// Assign FE to cells
hp::DoFHandler<2>::active_cell_iterator cell;
data_out.write_vtu (deallog.get_file_stream());
+ // Interpoalte solution
+ SolutionTransfer<2, Vector<double>, hp::DoFHandler<2> > solultion_trans(dof_handler);
+ solultion_trans.prepare_for_coarsening_and_refinement(solution);
+
+ triangulation.execute_coarsening_and_refinement ();
// Assign FE_Q_ to all cells
cell = dof_handler.begin_active();
for( ; cell != endc; ++cell){
cell->set_active_fe_index(0);
}
-
-
- // Interpoalte solution
- SolutionTransfer<2, Vector<double>, hp::DoFHandler<2> > solultion_trans(dof_handler);
- solultion_trans.prepare_for_coarsening_and_refinement(solution);
- triangulation.execute_coarsening_and_refinement ();
dof_handler.distribute_dofs (fe_collection);
+
Vector<double> new_solution(dof_handler.n_dofs());
solultion_trans.interpolate(solution, new_solution);