for (hp::DoFHandler<2>::active_cell_iterator cell =
dof_handler_patch.begin_active (); cell != dof_handler_patch.end ();
++cell, ++cell_index) {
dofs_per_cell = cell->get_fe ().dofs_per_cell;
cell_patch_dof_indices.resize (dofs_per_cell);
cell->get_dof_indices (cell_patch_dof_indices);
cell_dof_indices.resize (cell_list[cell_index]->get_fe
().dofs_per_cell);
cell_list[cell_index]->get_dof_indices (cell_dof_indices);
Are you sure that cell_list[cell_index] is the cell you are looking for?
I'm not entirely sure if Triangulation guarantees that the order of cells
is maintained when creating a mesh. In any case, it may be rotated, for
example, in which case the vertices of cell and cell_list[cell_index]
don't match any more.
I don't know whether that's the case, but it would be worth testing.
for (unsigned int i = 0; i < dofs_per_cell; ++i)
solution_patch_tmp (cell_patch_dof_indices[i]) = solution
(cell_dof_indices[i]);
Assuming the code works, this can probably be done in a simpler way by
using
cell_list[cell_index]->get_dof_values (solution, local_solution_values);
cell->set_dof_values (local_solution_values, solution_patch_tmp);
If you do this, you won't have to deal with DoF indices at all any more.
SolutionTransfer<2, Vector<double>, hp::DoFHandler<2> >
solution_transfer (dof_handler_patch);
for (hp::DoFHandler<2>::active_cell_iterator cell =
dof_handler_patch.begin_active (); cell != dof_handler_patch.end ();
++cell)
if (cell->material_id () == 1) {
cell->set_refine_flag ();
for (unsigned int face = 0; face <
GeometryInfo<2>::faces_per_cell; ++face)
if (!(cell->face (face)->at_boundary () || cell->face
(face)->has_children ()))
cell->neighbor (face)->set_refine_flag ();
}
triangulation_patch.prepare_coarsening_and_refinement ();
solution_transfer.prepare_for_pure_refinement ();
triangulation_patch.execute_coarsening_and_refinement ();
dof_handler_patch.distribute_dofs (fe_collection);
const unsigned int n_dofs_refined = dof_handler_patch.n_dofs ();
Vector<double> solution_patch (n_dofs_refined);
solution_transfer.refine_interpolate (solution_patch_tmp,
solution_patch);
constraints.clear ();
DoFTools::make_hanging_node_constraints (dof_handler_patch,
constraints); constraints.close ();
constraints.distribute (solution_patch);
For the extrapolation I can use SolutionTransfer, right? Unfortunately
the extrapolated solution looks very strange and does not correspond to
the solution on the old grid.
You should be able to. You could also first refine the mesh and then use
cell->set_dof_values_by_interpolation (...)
in the call above, which will write the result not to the current cell but
automatically interpolate it all the way down to its children.
As for what's going wrong: I would first check that indeed the cells that
you think match do indeed match, i.e. that in your loop cell and
cell_list[cell_index] refer to the same cell with the same vertices in the
same order.
Best
W.
