From: Wolfgang Bangerth Date: Tue, 14 Jul 2009 02:17:51 +0000 (+0000) Subject: Implement one way to fix up cells with distorted children. X-Git-Tag: v8.0.0~7484 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=f74be032eb0dd68cc6d1872640e7234bbfc655cc;p=dealii.git Implement one way to fix up cells with distorted children. git-svn-id: https://svn.dealii.org/trunk@19076 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/deal.II/include/grid/grid_tools.h b/deal.II/deal.II/include/grid/grid_tools.h index fd53c06208..d34e920e9c 100644 --- a/deal.II/deal.II/include/grid/grid_tools.h +++ b/deal.II/deal.II/include/grid/grid_tools.h @@ -2,7 +2,7 @@ // $Id$ // Version: $Name$ // -// Copyright (C) 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008 by the deal.II authors +// Copyright (C) 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009 by the deal.II authors // // This file is subject to QPL and may not be distributed // without copyright and license information. Please refer @@ -485,8 +485,9 @@ class GridTools */ template static - void partition_triangulation (const unsigned int n_partitions, - Triangulation &triangulation); + void + partition_triangulation (const unsigned int n_partitions, + Triangulation &triangulation); /** * This function does the same as the @@ -566,9 +567,10 @@ class GridTools */ template static - void partition_triangulation (const unsigned int n_partitions, - const SparsityPattern &cell_connection_graph, - Triangulation &triangulation); + void + partition_triangulation (const unsigned int n_partitions, + const SparsityPattern &cell_connection_graph, + Triangulation &triangulation); /** * For each active cell, return in the @@ -756,6 +758,37 @@ class GridTools create_union_triangulation (const Triangulation &triangulation_1, const Triangulation &triangulation_2, Triangulation &result); + + /** + * Given a triangulation and a + * list of cells whose children + * have become distorted as a + * result of mesh refinement, try + * to fix these cells up by + * moving the center node around. + * + * The function returns a list of + * cells with distorted children + * that couldn't be fixed up for + * whatever reason. The returned + * list is therefore a subset of + * the input argument. + * + * For a definition of the + * concept of distorted cells, + * see the + * @ref GlossDistorted "glossary entry". + * The first argument passed to the + * current function is typically + * the exception thrown by the + * Triangulation::execute_coarsening_and_refinement + * function. + */ + template + static + typename Triangulation::DistortedCellList + fix_up_distorted_child_cells (const typename Triangulation::DistortedCellList &distorted_cells, + Triangulation &triangulation); /** * Exception diff --git a/deal.II/deal.II/source/grid/grid_tools.cc b/deal.II/deal.II/source/grid/grid_tools.cc index a2b3e7f856..f93bfedac4 100644 --- a/deal.II/deal.II/source/grid/grid_tools.cc +++ b/deal.II/deal.II/source/grid/grid_tools.cc @@ -31,6 +31,8 @@ #include #include +#include + DEAL_II_NAMESPACE_OPEN @@ -1232,6 +1234,258 @@ GridTools::create_union_triangulation (const Triangulation &trian } +namespace internal +{ + namespace GridTools + { + namespace FixUpDistortedChildCells + { + template + double + objective_function (const typename dealii::Triangulation::cell_iterator &cell, + const Point &cell_mid_point) + { + // everything below is wrong + // if not for the following + // condition + Assert (cell->refinement_case() == RefinementCase::isotropic_refinement, + ExcNotImplemented()); + // first calculate the + // average jacobian + // determinant for the parent + // cell + Point parent_vertices + [GeometryInfo::vertices_per_cell]; + double parent_determinants + [GeometryInfo::vertices_per_cell]; + + for (unsigned int i=0; i::vertices_per_cell; ++i) + parent_vertices[i] = cell->vertex(i); + + GeometryInfo::jacobian_determinants_at_vertices (parent_vertices, + parent_determinants); + + const double average_parent_jacobian + = std::accumulate (&parent_determinants[0], + &parent_determinants[GeometryInfo::vertices_per_cell], + 0.); + + // now do the same + // computation for the + // children where we use the + // given location for the + // cell mid point instead of + // the one the triangulation + // currently reports + Point child_vertices + [GeometryInfo::max_children_per_cell] + [GeometryInfo::vertices_per_cell]; + double child_determinants + [GeometryInfo::max_children_per_cell] + [GeometryInfo::vertices_per_cell]; + + for (unsigned int c=0; cn_children(); ++c) + for (unsigned int i=0; i::vertices_per_cell; ++i) + child_vertices[c][i] = cell->child(c)->vertex(i); + + // replace mid-cell + // vertex. note that for + // child i, the mid-cell + // vertex happens to have the + // number + // max_children_per_cell-i + for (unsigned int c=0; cn_children(); ++c) + child_vertices[c][GeometryInfo::max_children_per_cell-c-1] + = cell_mid_point; + + for (unsigned int c=0; cn_children(); ++c) + GeometryInfo::jacobian_determinants_at_vertices (child_vertices[c], + child_determinants[c]); + + // on a uniformly refined + // hypercube cell, the child + // jacobians should all be + // smaller by a factor of + // 2^dim than the ones of the + // parent. as a consequence, + // we'll use the squared + // deviation from this ideal + // value as an objective + // function + double objective = 0; + for (unsigned int c=0; cn_children(); ++c) + for (unsigned int i=0; i::vertices_per_cell; ++i) + objective += std::pow (child_determinants[c][i] - + average_parent_jacobian/std::pow(2.,1.*dim), + 2); + + return objective; + } + + } + } +} + + + +template +typename Triangulation::DistortedCellList +GridTools:: +fix_up_distorted_child_cells (const typename Triangulation::DistortedCellList &distorted_cells, + Triangulation &/*triangulation*/) +{ + typename Triangulation::DistortedCellList unfixable_subset; + + // loop over all cells that we have + // to fix up + for (typename std::list::active_cell_iterator>::const_iterator + cell_ptr = distorted_cells.distorted_cells.begin(); + cell_ptr != distorted_cells.distorted_cells.end(); ++cell_ptr) + { + const typename Triangulation::cell_iterator + cell = *cell_ptr; + + // right now we can only deal + // with cells that have been + // refined isotropically + // because that is the only + // case where we have a cell + // mid-point that can be moved + // around without having to + // consider boundary + // information + Assert (cell->has_children(), ExcInternalError()); + Assert (cell->refinement_case() == RefinementCase::isotropic_refinement, + ExcNotImplemented()); + + // get the current location of + // the cell mid-vertex: + Point cell_mid_point + = cell->child(0)->vertex (GeometryInfo::max_children_per_cell-1); + + // now do a few steepest + // descent steps to reduce the + // objective function + unsigned int iteration = 0; + do + { + // choose a step length + // that is initially 1/10 + // of the child cells' + // diameter, and a sequence + // whose sum does not + // converge (to avoid + // premature termination of + // the iteration) + const double step_length = cell->diameter() / 10 / (iteration + 1); + + // compute the objective + // function and its derivative + const double val = internal::GridTools::FixUpDistortedChildCells:: + objective_function (cell, cell_mid_point); + + Tensor<1,dim> gradient; + for (unsigned int d=0; d h; + h[d] = step_length/2; + gradient[d] = (internal::GridTools::FixUpDistortedChildCells:: + objective_function (cell, + cell_mid_point + h) + - + internal::GridTools::FixUpDistortedChildCells:: + objective_function (cell, + cell_mid_point - h)) + / + step_length; + } + + // so we need to go in + // direction -gradient. the + // optimal value of the + // objective function is + // zero, so assuming that + // the model is quadratic + // we would have to go + // -2*val/||gradient|| in + // this direction, make + // sure we go at most + // step_length into this + // direction + cell_mid_point -= std::min(2*val / (gradient*gradient), + step_length / gradient.norm()) * + gradient; + + ++iteration; + } + while (iteration < 10); + + + // verify that the new location + // is indeed better than the + // one before in terms of the + // minimal jacobian determinant + double old_min_jacobian, new_min_jacobian; + + for (unsigned int test=0; test<2; ++test) + { + Point child_vertices + [GeometryInfo::max_children_per_cell] + [GeometryInfo::vertices_per_cell]; + double child_determinants + [GeometryInfo::max_children_per_cell] + [GeometryInfo::vertices_per_cell]; + + if (test == 1) + for (unsigned int c=0; cn_children(); ++c) + for (unsigned int i=0; i::vertices_per_cell; ++i) + child_vertices[c][i] = cell->child(c)->vertex(i); + + // replace mid-cell + // vertex. note that for + // child i, the mid-cell + // vertex happens to have the + // number + // max_children_per_cell-i + for (unsigned int c=0; cn_children(); ++c) + child_vertices[c][GeometryInfo::max_children_per_cell-c-1] + = cell_mid_point; + + for (unsigned int c=0; cn_children(); ++c) + GeometryInfo::jacobian_determinants_at_vertices (child_vertices[c], + child_determinants[c]); + + double min = child_determinants[0][0]; + for (unsigned int c=0; cn_children(); ++c) + for (unsigned int i=0; i::vertices_per_cell; ++i) + min = std::min (min, + child_determinants[c][i]); + + if (test == 0) + old_min_jacobian = min; + else + new_min_jacobian = min; + } + + // if new minimum jacobian is + // better than before, and if + // it is positive, then set the + // new mid point. otherwise + // return this cell as one of + // those that can't apparently + // be fixed + if ((new_min_jacobian > old_min_jacobian) + && + (new_min_jacobian > 0)) + cell->child(0)->vertex (GeometryInfo::max_children_per_cell-1) + = cell_mid_point; + else + unfixable_subset.distorted_cells.push_back (cell); + } + + return unfixable_subset; +} + // explicit instantiations @@ -1318,6 +1572,13 @@ GridTools::create_union_triangulation (const Triangulation &t const Triangulation &triangulation_2, Triangulation &result); +template +Triangulation::DistortedCellList +GridTools:: +fix_up_distorted_child_cells (const Triangulation::DistortedCellList &distorted_cells, + Triangulation &triangulation); + + #if deal_II_dimension != 3 @@ -1342,7 +1603,6 @@ void GridTools::scale (const double, Triangulation &); - #endif diff --git a/deal.II/doc/doxygen/headers/glossary.h b/deal.II/doc/doxygen/headers/glossary.h index ce3fd1b445..d0b2eb29b3 100644 --- a/deal.II/doc/doxygen/headers/glossary.h +++ b/deal.II/doc/doxygen/headers/glossary.h @@ -159,6 +159,11 @@ * Triangulation::execute_coarsening_and_refinement can then decide * what to do with this situation. * + * One way to deal with this problem is to use the + * GridTools::fix_up_distorted_child_cells function that attempts to + * fix up exactly these cells if possible by moving around the node at + * the center of the cell. + * * *
@anchor GlossFaceOrientation Face orientation
*
In a triangulation, the normal vector to a face diff --git a/tests/bits/distorted_cells_04.cc b/tests/bits/distorted_cells_04.cc new file mode 100644 index 0000000000..c82d443cb6 --- /dev/null +++ b/tests/bits/distorted_cells_04.cc @@ -0,0 +1,122 @@ +//---------------------------- distorted_cells_04.cc --------------------------- +// $Id$ +// Version: $Name$ +// +// Copyright (C) 2003, 2004, 2005, 2009 by the deal.II authors +// +// This file is subject to QPL and may not be distributed +// without copyright and license information. Please refer +// to the file deal.II/doc/license.html for the text and +// further information on this license. +// +//---------------------------- distorted_cells_04.cc --------------------------- + + +// like _03, but catch the exception and pass it to GridTools::fix_up_distorted_child_cells + +#include "../tests.h" +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include + + +template +class MyBoundary : public Boundary +{ + virtual Point + get_new_point_on_line (const typename Triangulation::line_iterator &line) const + { + deallog << "Finding point between " + << line->vertex(0) << " and " + << line->vertex(1) << std::endl; + + // in 2d, find a point that + // lies on the opposite side + // of the quad. in 3d, choose + // the midpoint of the edge + if (dim == 2) + return Point(0,0.75); + else + return (line->vertex(0) + line->vertex(1)) / 2; + } + + virtual Point + get_new_point_on_quad (const typename Triangulation::quad_iterator &quad) const + { + deallog << "Finding point between " + << quad->vertex(0) << " and " + << quad->vertex(1) << " and " + << quad->vertex(2) << " and " + << quad->vertex(3) << std::endl; + + return Point(0,0,.75); + } +}; + + + +template +void check () +{ + MyBoundary my_boundary; + + // create a single square/cube + Triangulation coarse_grid; + GridGenerator::hyper_cube (coarse_grid, -1, 1); + + // set bottom face to use MyBoundary + for (unsigned int f=0; f::faces_per_cell; ++f) + if (coarse_grid.begin_active()->face(f)->center()[dim-1] == -1) + coarse_grid.begin_active()->face(f)->set_boundary_indicator (1); + coarse_grid.set_boundary (1, my_boundary); + + // now try to refine this one + // cell. we should get an exception + try + { + coarse_grid.begin_active()->set_refine_flag (); + coarse_grid.execute_coarsening_and_refinement (); + } + catch (typename Triangulation::DistortedCellList &dcv) + { + typename Triangulation::DistortedCellList + subset = GridTools::fix_up_distorted_child_cells (dcv, + coarse_grid); + Assert (subset.distorted_cells.size() == 0, + ExcInternalError()); + } + + Assert (coarse_grid.n_levels() == 2, ExcInternalError()); + Assert (coarse_grid.n_active_cells() == 1< (); + check<3> (); +} + + + diff --git a/tests/bits/distorted_cells_04/2d b/tests/bits/distorted_cells_04/2d new file mode 100644 index 0000000000..1a53889d21 --- /dev/null +++ b/tests/bits/distorted_cells_04/2d @@ -0,0 +1,6 @@ +4 1 0 0 0 +1 0 0 0 +2 1 0 0 +3 0 1 0 +4 1 1 0 +1 0 quad 1 2 3 4 diff --git a/tests/bits/distorted_cells_04/3d b/tests/bits/distorted_cells_04/3d new file mode 100644 index 0000000000..dd37099c80 --- /dev/null +++ b/tests/bits/distorted_cells_04/3d @@ -0,0 +1,10 @@ +8 1 0 0 0 +1 0 0 0 +2 1 0 0 +3 0 1 0 +4 1 1 0 +5 0 0 1 +6 1 0 1 +7 0 1 1 +8 1 1 1 +1 0 hex 1 2 3 4 5 6 7 8 diff --git a/tests/bits/distorted_cells_04/cmp/generic b/tests/bits/distorted_cells_04/cmp/generic new file mode 100644 index 0000000000..bba3339d5d --- /dev/null +++ b/tests/bits/distorted_cells_04/cmp/generic @@ -0,0 +1,223 @@ + +DEAL::Finding point between -1.00000 -1.00000 and 1.00000 -1.00000 +-1.00000 -1.00000 1 0 +0.00000 0.750000 1 0 +-3.50619e-10 0.846953 1 0 +-1.00000 0.00000 1 0 +-1.00000 -1.00000 1 0 + + +0.00000 0.750000 1 0 +1.00000 -1.00000 1 0 +1.00000 0.00000 1 0 +-3.50619e-10 0.846953 1 0 +0.00000 0.750000 1 0 + + +-1.00000 0.00000 1 0 +-3.50619e-10 0.846953 1 0 +0.00000 1.00000 1 0 +-1.00000 1.00000 1 0 +-1.00000 0.00000 1 0 + + +-3.50619e-10 0.846953 1 0 +1.00000 0.00000 1 0 +1.00000 1.00000 1 0 +0.00000 1.00000 1 0 +-3.50619e-10 0.846953 1 0 + + +DEAL::Finding point between -1.00000 -1.00000 -1.00000 and 1.00000 -1.00000 -1.00000 and -1.00000 1.00000 -1.00000 and 1.00000 1.00000 -1.00000 +-1.00000 -1.00000 -1.00000 1 0 +0.00000 -1.00000 -1.00000 1 0 +0.00000 -1.00000 0.00000 1 0 +-1.00000 -1.00000 0.00000 1 0 +-1.00000 -1.00000 -1.00000 1 0 + +-1.00000 0.00000 -1.00000 1 0 +0.00000 0.00000 0.750000 1 0 +1.07919e-09 -5.43281e-10 0.876731 1 0 +-1.00000 0.00000 0.00000 1 0 +-1.00000 0.00000 -1.00000 1 0 + +-1.00000 -1.00000 -1.00000 1 0 +-1.00000 0.00000 -1.00000 1 0 + +0.00000 -1.00000 -1.00000 1 0 +0.00000 0.00000 0.750000 1 0 + +0.00000 -1.00000 0.00000 1 0 +1.07919e-09 -5.43281e-10 0.876731 1 0 + +-1.00000 -1.00000 0.00000 1 0 +-1.00000 0.00000 0.00000 1 0 + +0.00000 -1.00000 -1.00000 1 0 +1.00000 -1.00000 -1.00000 1 0 +1.00000 -1.00000 0.00000 1 0 +0.00000 -1.00000 0.00000 1 0 +0.00000 -1.00000 -1.00000 1 0 + +0.00000 0.00000 0.750000 1 0 +1.00000 0.00000 -1.00000 1 0 +1.00000 0.00000 0.00000 1 0 +1.07919e-09 -5.43281e-10 0.876731 1 0 +0.00000 0.00000 0.750000 1 0 + +0.00000 -1.00000 -1.00000 1 0 +0.00000 0.00000 0.750000 1 0 + +1.00000 -1.00000 -1.00000 1 0 +1.00000 0.00000 -1.00000 1 0 + +1.00000 -1.00000 0.00000 1 0 +1.00000 0.00000 0.00000 1 0 + +0.00000 -1.00000 0.00000 1 0 +1.07919e-09 -5.43281e-10 0.876731 1 0 + +-1.00000 0.00000 -1.00000 1 0 +0.00000 0.00000 0.750000 1 0 +1.07919e-09 -5.43281e-10 0.876731 1 0 +-1.00000 0.00000 0.00000 1 0 +-1.00000 0.00000 -1.00000 1 0 + +-1.00000 1.00000 -1.00000 1 0 +0.00000 1.00000 -1.00000 1 0 +0.00000 1.00000 0.00000 1 0 +-1.00000 1.00000 0.00000 1 0 +-1.00000 1.00000 -1.00000 1 0 + +-1.00000 0.00000 -1.00000 1 0 +-1.00000 1.00000 -1.00000 1 0 + +0.00000 0.00000 0.750000 1 0 +0.00000 1.00000 -1.00000 1 0 + +1.07919e-09 -5.43281e-10 0.876731 1 0 +0.00000 1.00000 0.00000 1 0 + +-1.00000 0.00000 0.00000 1 0 +-1.00000 1.00000 0.00000 1 0 + +0.00000 0.00000 0.750000 1 0 +1.00000 0.00000 -1.00000 1 0 +1.00000 0.00000 0.00000 1 0 +1.07919e-09 -5.43281e-10 0.876731 1 0 +0.00000 0.00000 0.750000 1 0 + +0.00000 1.00000 -1.00000 1 0 +1.00000 1.00000 -1.00000 1 0 +1.00000 1.00000 0.00000 1 0 +0.00000 1.00000 0.00000 1 0 +0.00000 1.00000 -1.00000 1 0 + +0.00000 0.00000 0.750000 1 0 +0.00000 1.00000 -1.00000 1 0 + +1.00000 0.00000 -1.00000 1 0 +1.00000 1.00000 -1.00000 1 0 + +1.00000 0.00000 0.00000 1 0 +1.00000 1.00000 0.00000 1 0 + +1.07919e-09 -5.43281e-10 0.876731 1 0 +0.00000 1.00000 0.00000 1 0 + +-1.00000 -1.00000 0.00000 1 0 +0.00000 -1.00000 0.00000 1 0 +0.00000 -1.00000 1.00000 1 0 +-1.00000 -1.00000 1.00000 1 0 +-1.00000 -1.00000 0.00000 1 0 + +-1.00000 0.00000 0.00000 1 0 +1.07919e-09 -5.43281e-10 0.876731 1 0 +0.00000 0.00000 1.00000 1 0 +-1.00000 0.00000 1.00000 1 0 +-1.00000 0.00000 0.00000 1 0 + +-1.00000 -1.00000 0.00000 1 0 +-1.00000 0.00000 0.00000 1 0 + +0.00000 -1.00000 0.00000 1 0 +1.07919e-09 -5.43281e-10 0.876731 1 0 + +0.00000 -1.00000 1.00000 1 0 +0.00000 0.00000 1.00000 1 0 + +-1.00000 -1.00000 1.00000 1 0 +-1.00000 0.00000 1.00000 1 0 + +0.00000 -1.00000 0.00000 1 0 +1.00000 -1.00000 0.00000 1 0 +1.00000 -1.00000 1.00000 1 0 +0.00000 -1.00000 1.00000 1 0 +0.00000 -1.00000 0.00000 1 0 + +1.07919e-09 -5.43281e-10 0.876731 1 0 +1.00000 0.00000 0.00000 1 0 +1.00000 0.00000 1.00000 1 0 +0.00000 0.00000 1.00000 1 0 +1.07919e-09 -5.43281e-10 0.876731 1 0 + +0.00000 -1.00000 0.00000 1 0 +1.07919e-09 -5.43281e-10 0.876731 1 0 + +1.00000 -1.00000 0.00000 1 0 +1.00000 0.00000 0.00000 1 0 + +1.00000 -1.00000 1.00000 1 0 +1.00000 0.00000 1.00000 1 0 + +0.00000 -1.00000 1.00000 1 0 +0.00000 0.00000 1.00000 1 0 + +-1.00000 0.00000 0.00000 1 0 +1.07919e-09 -5.43281e-10 0.876731 1 0 +0.00000 0.00000 1.00000 1 0 +-1.00000 0.00000 1.00000 1 0 +-1.00000 0.00000 0.00000 1 0 + +-1.00000 1.00000 0.00000 1 0 +0.00000 1.00000 0.00000 1 0 +0.00000 1.00000 1.00000 1 0 +-1.00000 1.00000 1.00000 1 0 +-1.00000 1.00000 0.00000 1 0 + +-1.00000 0.00000 0.00000 1 0 +-1.00000 1.00000 0.00000 1 0 + +1.07919e-09 -5.43281e-10 0.876731 1 0 +0.00000 1.00000 0.00000 1 0 + +0.00000 0.00000 1.00000 1 0 +0.00000 1.00000 1.00000 1 0 + +-1.00000 0.00000 1.00000 1 0 +-1.00000 1.00000 1.00000 1 0 + +1.07919e-09 -5.43281e-10 0.876731 1 0 +1.00000 0.00000 0.00000 1 0 +1.00000 0.00000 1.00000 1 0 +0.00000 0.00000 1.00000 1 0 +1.07919e-09 -5.43281e-10 0.876731 1 0 + +0.00000 1.00000 0.00000 1 0 +1.00000 1.00000 0.00000 1 0 +1.00000 1.00000 1.00000 1 0 +0.00000 1.00000 1.00000 1 0 +0.00000 1.00000 0.00000 1 0 + +1.07919e-09 -5.43281e-10 0.876731 1 0 +0.00000 1.00000 0.00000 1 0 + +1.00000 0.00000 0.00000 1 0 +1.00000 1.00000 0.00000 1 0 + +1.00000 0.00000 1.00000 1 0 +1.00000 1.00000 1.00000 1 0 + +0.00000 0.00000 1.00000 1 0 +0.00000 1.00000 1.00000 1 0 +