From 180e842e48cdf38602a8fcf1989d85cdb0e049a7 Mon Sep 17 00:00:00 2001 From: heltai Date: Tue, 28 Jan 2014 14:42:29 +0000 Subject: [PATCH] Fixed two distorted cells tests. git-svn-id: https://svn.dealii.org/branches/branch_manifold_id@32333 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/include/deal.II/grid/manifold.h | 19 +++++ deal.II/source/grid/grid_tools.cc | 12 +++- deal.II/source/grid/manifold.cc | 9 +++ tests/bits/distorted_cells_03.cc | 39 +--------- tests/bits/distorted_cells_04.cc | 73 +------------------ tests/bits/my_boundary.h | 95 +++++++++++++++++++++++++ 6 files changed, 134 insertions(+), 113 deletions(-) create mode 100644 tests/bits/my_boundary.h diff --git a/deal.II/include/deal.II/grid/manifold.h b/deal.II/include/deal.II/grid/manifold.h index 955962b69c..48016a8286 100644 --- a/deal.II/include/deal.II/grid/manifold.h +++ b/deal.II/include/deal.II/grid/manifold.h @@ -83,6 +83,25 @@ public: Point get_new_point(const std::vector > &surrounding_points, const std::vector &weights) const = 0; + /** + * Given a point which lies close to the given manifold, it modifies + * it and projects it to manifold itself. This function is used in + * some mesh smoothing algorithms that try to move around points in + * order to improve the mesh quality but need to ensure that points + * that were on the manfold remain on the manifold. + * + * Derived classes do not need to implement this function unless + * mesh smoothing algorithms are used with a particular boundary + * object. The default implementation of this function calls the + * method with the same name, but without surrounding_points. + * + * This version of the projection function can be overwritten when + * only an approximation of the manifold is available, and this + * approximation requires surrounding points to be defined. + */ + virtual + Point project_to_manifold (const std::vector > &surrounding_points, + const Point &candidate) const; /** * Given a point which lies close to the given manifold, it modifies diff --git a/deal.II/source/grid/grid_tools.cc b/deal.II/source/grid/grid_tools.cc index e83cc14f5f..b385bb51b2 100644 --- a/deal.II/source/grid/grid_tools.cc +++ b/deal.II/source/grid/grid_tools.cc @@ -1751,6 +1751,9 @@ next_cell: const unsigned int structdim = Iterator::AccessorType::structure_dimension; const unsigned int spacedim = Iterator::AccessorType::space_dimension; + std::vector > vertices(GeometryInfo::vertices_per_cell); + for(unsigned int i=0; ivertex(i); // right now we can only deal // with cells that have been @@ -1816,10 +1819,12 @@ next_cell: else gradient[d] = ((objective_function (object, - manifold->project_to_manifold(object_mid_point + h)) + manifold->project_to_manifold(vertices, + object_mid_point + h)) - objective_function (object, - manifold->project_to_manifold(object_mid_point - h))) + manifold->project_to_manifold(vertices, + object_mid_point - h))) / eps); } @@ -1853,7 +1858,8 @@ next_cell: gradient; if (respect_manifold == true) - object_mid_point = manifold->project_to_manifold(object_mid_point); + object_mid_point = manifold->project_to_manifold(vertices, + object_mid_point); // compute current value of the // objective function diff --git a/deal.II/source/grid/manifold.cc b/deal.II/source/grid/manifold.cc index 60f6fb8f25..325430448f 100644 --- a/deal.II/source/grid/manifold.cc +++ b/deal.II/source/grid/manifold.cc @@ -40,6 +40,15 @@ Manifold::project_to_manifold (const Point &candidate) const return candidate; } +template +Point +Manifold::project_to_manifold (const std::vector > &, + const Point &candidate) const +{ + return project_to_manifold(candidate); +} + + template void Manifold::get_normals_at_points(const std::vector > &points, diff --git a/tests/bits/distorted_cells_03.cc b/tests/bits/distorted_cells_03.cc index 4d32febe71..4ab1d099e2 100644 --- a/tests/bits/distorted_cells_03.cc +++ b/tests/bits/distorted_cells_03.cc @@ -23,6 +23,7 @@ // distorted cells resulting from refinement #include "../tests.h" +#include "my_boundary.h" #include #include #include @@ -38,44 +39,6 @@ #include - -template -class MyBoundary : public Manifold -{ - virtual Point - get_new_point (const std::vector > &v, - const std::vector &ws) const - { - switch(v.size()) { - case 2: - deallog << "Finding point between " - << v[0] << " and " - << v[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 (v[0] + v[1]) / 2; - break; - default: - deallog << "Finding point between " - << v[0] << " and " - << v[1] << " and " - << v[2] << " and " - << v[3] << std::endl; - - return Point(0,0,.75); - break; - } - } -}; - - - template void check () { diff --git a/tests/bits/distorted_cells_04.cc b/tests/bits/distorted_cells_04.cc index 715cbd5759..1810e3cde2 100644 --- a/tests/bits/distorted_cells_04.cc +++ b/tests/bits/distorted_cells_04.cc @@ -34,78 +34,7 @@ #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); - } - - virtual - Point - project_to_surface (const typename Triangulation::line_iterator &, - const Point &p) const - { - deallog << "Projecting line point " << p << std::endl; - - if (dim == 2) - return Point(p[0], 0.75-1.75*std::fabs(p[0])); - else - return p; - } - - virtual - Point - project_to_surface (const typename Triangulation::quad_iterator &, - const Point &p) const - { - deallog << "Projecting quad point " << p << std::endl; - - Assert (dim == 3, ExcInternalError()); - - return Point(p[0], p[1], - 0.75-1.75*std::max(std::fabs(p[0]), - std::fabs(p[1]))); - } - - virtual - Point - project_to_surface (const typename Triangulation::hex_iterator &, - const Point &) const - { - Assert (false, ExcInternalError()); - return Point(); - } -}; - - +#include "my_boundary.h" template void check () diff --git a/tests/bits/my_boundary.h b/tests/bits/my_boundary.h new file mode 100644 index 0000000000..a251d30406 --- /dev/null +++ b/tests/bits/my_boundary.h @@ -0,0 +1,95 @@ +// --------------------------------------------------------------------- +// $Id: distorted_cells_04.cc 32209 2014-01-15 00:25:27Z heltai $ +// +// Copyright (C) 2003 - 2013 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + + + +// like _03, but catch the exception and pass it to GridTools::fix_up_distorted_child_cells + +#include "../tests.h" +#include + +template +class MyBoundary : public Manifold +{ + virtual Point + get_new_point (const std::vector > &v, + const std::vector &ws) const + { + switch(v.size()) { + case 2: + deallog << "Finding point between " + << v[0] << " and " + << v[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 (v[0] + v[1]) / 2; + break; + default: + deallog << "Finding point between " + << v[0] << " and " + << v[1] << " and " + << v[2] << " and " + << v[3] << std::endl; + + return Point(0,0,.75); + break; + } + } + + virtual + Point + project_to_manifold (const std::vector > &v, + const Point &p) const + { + switch(v.size()) { + case 2: + { + deallog << "Projecting line point " << p << std::endl; + if (dim == 2) + return Point(p[0], 0.75-1.75*std::fabs(p[0])); + else + return p; + } + break; + case 4: + { + deallog << "Projecting quad point " << p << std::endl; + + Assert (dim == 3, ExcInternalError()); + + return Point(p[0], p[1], + 0.75-1.75*std::max(std::fabs(p[0]), + std::fabs(p[1]))); + } + break; + default: + { + Assert (false, ExcInternalError()); + return Point(); + } + } + } +}; + + + + -- 2.39.5