From d4cb9852e2adaad7435072a3f7f0e58e48a75415 Mon Sep 17 00:00:00 2001 From: heister Date: Tue, 4 Feb 2014 18:43:20 +0000 Subject: [PATCH] make GridTools::transform work with hanging nodes git-svn-id: https://svn.dealii.org/trunk@32402 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/doc/news/changes.h | 5 ++ deal.II/include/deal.II/grid/grid_tools.h | 67 +++++++++++++++++------ 2 files changed, 56 insertions(+), 16 deletions(-) diff --git a/deal.II/doc/news/changes.h b/deal.II/doc/news/changes.h index be608b1720..6aba67d558 100644 --- a/deal.II/doc/news/changes.h +++ b/deal.II/doc/news/changes.h @@ -124,6 +124,11 @@ inconvenience this causes.

Specific improvements

    +
  1. Changed: GridTools::transform() can now deal with meshes with hanging nodes. +
    + (Timo Heister, 2014/02/04) +
  2. +
  3. Fixed: Calling FEValuesViews::Vector::get_function_curls() computed wrong results in some cases (see https://code.google.com/p/dealii/issues/detail?id=182). This is now fixed. diff --git a/deal.II/include/deal.II/grid/grid_tools.h b/deal.II/include/deal.II/grid/grid_tools.h index 5f43b3af9e..e8de5d95ab 100644 --- a/deal.II/include/deal.II/grid/grid_tools.h +++ b/deal.II/include/deal.II/grid/grid_tools.h @@ -171,15 +171,7 @@ namespace GridTools /** * Transform the vertices of the given * triangulation by applying the - * function object provided as first argument to all its vertices. Since - * the internal consistency of a - * triangulation can only be guaranteed - * if the transformation is applied to - * the vertices of only one level of - * hierarchically refined cells, this - * function may only be used if all cells - * of the triangulation are on the same - * refinement level. + * function object provided as first argument to all its vertices. * * The transformation given as * argument is used to transform @@ -194,6 +186,13 @@ namespace GridTools * return value have to be of * type Point. * + * Note: if you are using a parallel::distributed::Triangulation you will have + * hanging nodes in your local Triangulation even if your "global" mesh has + * no hanging nodes. This will cause issues with wrong positioning of hanging + * nodes in ghost cells. The active cells will be correct, but keep in mind that + * computations like KellyErrorEstimator will give wrong answers. A safe bet is + * to use this function prior to any refinement in parallel. + * * This function is used in the * "Possibilities for extensions" section * of step-38. It is also used in step-49. @@ -1204,13 +1203,6 @@ namespace GridTools void transform (const Predicate &predicate, Triangulation &triangulation) { - // ensure that all the cells of the - // triangulation are on the same level - Assert (triangulation.n_levels() == - static_cast(triangulation.begin_active()->level()+1), - ExcMessage ("Not all cells of this triangulation are at the same " - "refinement level, as is required for this function.")); - std::vector treated_vertices (triangulation.n_vertices(), false); @@ -1232,6 +1224,49 @@ namespace GridTools // and mark it as treated treated_vertices[cell->vertex_index(v)] = true; }; + + + // now fix any vertices on hanging nodes so that we don't create any holes + if (dim==2) + { + typename Triangulation::active_cell_iterator + cell = triangulation.begin_active(), + endc = triangulation.end(); + for (; cell!=endc; ++cell) + for (unsigned int face=0; face::faces_per_cell; ++face) + if (cell->face(face)->has_children() && + !cell->face(face)->at_boundary()) + { // this line has children + cell->face(face)->child(0)->vertex(1) + = (cell->face(face)->vertex(0) + + cell->face(face)->vertex(1)) / 2; + } + } + else if (dim==3) + { + typename Triangulation::active_cell_iterator + cell = triangulation.begin_active(), + endc = triangulation.end(); + for (; cell!=endc; ++cell) + for (unsigned int face=0; face::faces_per_cell; ++face) + if (cell->face(face)->has_children() && + !cell->face(face)->at_boundary()) + { // this face has hanging nodes + cell->face(face)->child(0)->vertex(1) + = (cell->face(face)->vertex(0) + cell->face(face)->vertex(1)) / 2.0; + cell->face(face)->child(0)->vertex(2) + = (cell->face(face)->vertex(0) + cell->face(face)->vertex(2)) / 2.0; + cell->face(face)->child(1)->vertex(3) + = (cell->face(face)->vertex(1) + cell->face(face)->vertex(3)) / 2.0; + cell->face(face)->child(2)->vertex(3) + = (cell->face(face)->vertex(2) + cell->face(face)->vertex(3)) / 2.0; + + // center of the face + cell->face(face)->child(0)->vertex(3) + = (cell->face(face)->vertex(0) + cell->face(face)->vertex(1) + + cell->face(face)->vertex(2) + cell->face(face)->vertex(3)) / 4.0; + } + } } -- 2.39.5