From: bangerth Date: Wed, 30 May 2012 21:11:51 +0000 (+0000) Subject: Improve an algorithm in MappingQ::transform_real_to_unit_cell. X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=b39f284f04307e4b723f482adf0c4f35a4b4e423;p=dealii-svn.git Improve an algorithm in MappingQ::transform_real_to_unit_cell. git-svn-id: https://svn.dealii.org/trunk@25578 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/doc/news/changes.h b/deal.II/doc/news/changes.h index be820c9c3b..8a4f303a51 100644 --- a/deal.II/doc/news/changes.h +++ b/deal.II/doc/news/changes.h @@ -201,13 +201,18 @@ enabled due to a missing include file in file

Specific improvements

    +
  1. Fixed: A problem in MappingQ::transform_real_to_unit_cell +that sometimes led the algorithm in this function to abort. +
    +(Wolfgang Bangerth 2012/05/30) +
  2. New: The function DataOutInterface::write_pvd_record can be used to provide Paraview with metadata that describes which time in a simulation a particular output file corresponds to.
    (Marco Engelhard 2012/05/30) -
  3. Fixed: Bug in 3d with hanging nodes in GridTools::find_cells_adjacent_to_vertex() +
  4. Fixed: A bug in 3d with hanging nodes in GridTools::find_cells_adjacent_to_vertex() that caused find_active_cell_around_point() to fail in those cases.
    (Timo Heister, Wolfgang Bangerth 2012/05/30) diff --git a/deal.II/source/fe/mapping_q.cc b/deal.II/source/fe/mapping_q.cc index a6e9222bf2..0beb8330d2 100644 --- a/deal.II/source/fe/mapping_q.cc +++ b/deal.II/source/fe/mapping_q.cc @@ -1445,9 +1445,10 @@ MappingQ:: transform_real_to_unit_cell (const typename Triangulation::cell_iterator &cell, const Point &p) const { - // first a Newton iteration based - // on a Q1 mapping - Point p_unit = MappingQ1::transform_real_to_unit_cell(cell, p); + // first a Newton iteration based on a Q1 + // mapping. + Point p_unit = + MappingQ1::transform_real_to_unit_cell(cell, p); // then a Newton iteration based on the // full MappingQ if we need this. note that @@ -1459,6 +1460,19 @@ transform_real_to_unit_cell (const typename Triangulation::cell_it use_mapping_q_on_all_cells || (dim!=spacedim) ) { + // use the full mapping. in case the + // function above should have given us + // something back that lies outside the + // unit cell (that might happen because + // we may have given a point 'p' that + // lies inside the cell with the higher + // order mapping, but outside the + // Q1-mapped reference cell), then + // project it back into the reference + // cell in hopes that this gives a + // better starting point to the + // following iteration + p_unit = GeometryInfo::project_to_unit_cell(p_unit); const Quadrature point_quadrature(p_unit); diff --git a/deal.II/source/fe/mapping_q1.cc b/deal.II/source/fe/mapping_q1.cc index b3d683b0b5..87485d368c 100644 --- a/deal.II/source/fe/mapping_q1.cc +++ b/deal.II/source/fe/mapping_q1.cc @@ -1933,10 +1933,8 @@ transform_real_to_unit_cell_internal_codim1 // eps. Only if the first condition // failed, loop will have been // increased and tested, and thus - // havereached the limit. - AssertThrow(loop(); -// test_real_to_unit_cell<2,3>(); return 0; } diff --git a/tests/fe/mapping_real_to_unit_q4_curved_codim.cc b/tests/fe/mapping_real_to_unit_q4_curved_codim.cc new file mode 100644 index 0000000000..2349ba77e1 --- /dev/null +++ b/tests/fe/mapping_real_to_unit_q4_curved_codim.cc @@ -0,0 +1,121 @@ +//----------------------------------------------------------------------------- +// $Id: data_out_base_pvd.cc 25569 2012-05-30 12:53:31Z bangerth $ +// Version: $Name$ +// +// Copyright (C) 2006, 2007, 2010, 2012 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. +// +//----------------------------------------------------------------------------- + + +// on a somewhat deformed cube, verify that if we push forward a bunch +// of points from the reference to the real cell and then call +// Mapping::transform_unit_to_real_cell that we get the same point as +// we had in the beginning. +// +// like in the _q4_straight test, we use a Q4 mapping but this time we +// actually curve one boundary of the cell which ensures that the +// mapping is really higher order than just Q1 + +#include "../tests.h" + +#include +#include +#include +#include +#include + + +template +void test_real_to_unit_cell() +{ + deallog << "dim=" << dim << ", spacedim=" << spacedim << std::endl; + + // define a boundary that fits the + // the vertices of the hyper cube + // we're going to create below + HyperBallBoundary boundary (Point(), + std::sqrt(1.*dim)); + + Triangulation triangulation; + GridGenerator::hyper_cube (triangulation, -1, 1); + + // set the boundary indicator for + // one face of the single cell + triangulation.set_boundary (1, boundary); + triangulation.begin_active()->face(0)->set_boundary_indicator (1); + + const unsigned int n_points = 5; + std::vector< Point > unit_points(Utilities::fixed_power(n_points)); + + switch (dim) + { + case 1: + for (unsigned int x=0; x map(4); + + // work with this cell (unlike the + // _q1 test where we move vertices) + typename Triangulation::active_cell_iterator + cell = triangulation.begin_active(); + for (unsigned int i=0; i p = map.transform_unit_to_real_cell(cell,unit_points[i]); + const Point p_unit = map.transform_real_to_unit_cell(cell,p); + Assert (unit_points[i].distance(p_unit) < 1e-10, + ExcInternalError()); + } + deallog << "OK" << std::endl; +} + + +int +main() +{ + std::ofstream logfile ("mapping_real_to_unit_q4_curved_codim/output"); + deallog.attach(logfile); + deallog.depth_console(0); + deallog.threshold_double(1.e-10); + + + test_real_to_unit_cell<2,3>(); + + return 0; +} + + + diff --git a/tests/fe/mapping_real_to_unit_q4_curved_codim/cmp/generic b/tests/fe/mapping_real_to_unit_q4_curved_codim/cmp/generic new file mode 100644 index 0000000000..4e66a74d6e --- /dev/null +++ b/tests/fe/mapping_real_to_unit_q4_curved_codim/cmp/generic @@ -0,0 +1,3 @@ + +DEAL::dim=2, spacedim=3 +DEAL::OK