Hello,
I am extracting a part of the triangulation and the corresponding
solution to a new grid, where cell_list containes the cells to be
extracted. Then I want to refine some cells and extrapolate the my
previously computed solution onto the new grid. This is mainly done by
the following code:
const unsigned int n_cells = cell_list.size ();
FullMatrix<double> cell_matrix (n_cells,
GeometryInfo<2>::vertices_per_cell);
for (unsigned int i = 0; i < n_cells; ++i)
for (unsigned int j = 0; j < GeometryInfo<2>::vertices_per_cell; ++j)
cell_matrix (i, j) = -1;
std::vector<Point<2> > vertices;
vertices.reserve (20);
Point<2> vertex;
unsigned int cell_index = 0;
unsigned int n_vertices;
for
(std::vector<hp::DoFHandler<2>::active_cell_iterator>::const_iterator
cell_ptr = cell_list.begin (); cell_ptr != cell_list.end (); ++cell_ptr,
++cell_index)
for (unsigned int i = 0; i < GeometryInfo<2>::vertices_per_cell;
++i) {
n_vertices = vertices.size ();
vertex = (*cell_ptr)->vertex (i);
for (unsigned int j = 0; j < n_vertices; ++j)
if (vertex == vertices[j]) {
cell_matrix (cell_index, i) = j;
break;
}
if (cell_matrix (cell_index, i) == -1) {
cell_matrix (cell_index, i) = n_vertices;
vertices.push_back (vertex);
}
}
std::vector<CellData<2> > cells (n_cells, CellData<2> ());
for (unsigned int i = 0; i < GeometryInfo<2>::vertices_per_cell; ++i)
cells[0].vertices[i] = cell_matrix (0, i);
cells[0].material_id = 1;
for (unsigned int i = 1; i < n_cells; ++i) {
for (unsigned int j = 0; j < GeometryInfo<2>::vertices_per_cell; ++j)
cells[i].vertices[j] = cell_matrix (i, j);
cells[i].material_id = 0;
}
Triangulation<2> triangulation_patch;
triangulation_patch.create_triangulation (vertices, cells, SubCellData
());
cell_index = 0;
bool pure_neumann_boundary = true;
hp::DoFHandler<2> dof_handler_patch (triangulation_patch);
for (hp::DoFHandler<2>::active_cell_iterator cell =
dof_handler_patch.begin_active (); cell != dof_handler_patch.end ();
++cell, ++cell_index) {
cell->set_active_fe_index (cell_list[cell_index]->active_fe_index ());
if (cell_list[cell_index]->at_boundary ())
for (unsigned int face = 0; face <
GeometryInfo<2>::faces_per_cell; ++face)
if (cell_list[cell_index]->face (face)->at_boundary ()) {
cell->face (face)->set_boundary_indicator (1);
pure_neumann_boundary = false;
}
}
dof_handler_patch.distribute_dofs (fe_collection);
cell_index = 0;
const unsigned int n_dofs = dof_handler_patch.n_dofs ();
std::vector<unsigned int> cell_dof_indices;
std::vector<unsigned int> cell_patch_dof_indices;
unsigned int dofs_per_cell;
Vector<double> solution_patch_tmp (n_dofs);
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);
for (unsigned int i = 0; i < dofs_per_cell; ++i)
solution_patch_tmp (cell_patch_dof_indices[i]) = solution
(cell_dof_indices[i]);
}
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.
What am I doing wrong here?
Thanks,
Markus
_______________________________________________
dealii mailing list http://poisson.dealii.org/mailman/listinfo/dealii
