Skip to content

Commit

Permalink
Merge pull request #13278 from kronbichler/fix_real_to_unit
Browse files Browse the repository at this point in the history
Fix bug in initial condition for transform_real_to_unit_cell
  • Loading branch information
kronbichler committed Jan 22, 2022
2 parents 6bbe680 + 91b7121 commit 4e8af90
Show file tree
Hide file tree
Showing 6 changed files with 552 additions and 402 deletions.
5 changes: 5 additions & 0 deletions doc/news/changes/minor/20220121MartinKronbichler
Original file line number Diff line number Diff line change
@@ -0,0 +1,5 @@
Fixed: MappingQ::transform_points_real_to_unit_cell has been made more robust,
correcting an unsuitable initial guess for the Newton iteration to compute the
inverse transformation in the case of bilinear/trilinear cells.
<br>
(Bruno Blais, Martin Kronbichler, 2022/01/22)
64 changes: 49 additions & 15 deletions include/deal.II/fe/mapping_q_internal.h
Original file line number Diff line number Diff line change
Expand Up @@ -482,6 +482,10 @@ namespace internal
const std::vector<unsigned int> & renumber,
const bool print_iterations_to_deallog = false)
{
if (print_iterations_to_deallog)
deallog << "Start MappingQ::do_transform_real_to_unit_cell for real "
<< "point [ " << p << " ] " << std::endl;

AssertDimension(points.size(),
Utilities::pow(polynomials_1d.size(), dim));

Expand Down Expand Up @@ -548,8 +552,8 @@ namespace internal
const unsigned int newton_iteration_limit = 20;

Point<dim, Number> invalid_point;
invalid_point[0] = std::numeric_limits<double>::infinity();
bool try_project_to_unit_cell = false;
invalid_point[0] = std::numeric_limits<double>::infinity();
bool tried_project_to_unit_cell = false;

unsigned int newton_iteration = 0;
Number f_weighted_norm_square = 1.;
Expand Down Expand Up @@ -578,7 +582,7 @@ namespace internal
// back to the unit cell and go on with the Newton iteration
// from there. Since the outside case is unlikely, we can
// afford spending the extra effort at this place.
if (try_project_to_unit_cell == false)
if (tried_project_to_unit_cell == false)
{
p_unit = GeometryInfo<dim>::project_to_unit_cell(p_unit);
p_real = internal::evaluate_tensor_product_value_and_gradient(
Expand All @@ -590,7 +594,7 @@ namespace internal
f = p_real.first - p;
f_weighted_norm_square = 1.;
last_f_weighted_norm_square = 1;
try_project_to_unit_cell = true;
tried_project_to_unit_cell = true;
continue;
}
else
Expand All @@ -606,7 +610,7 @@ namespace internal
deallog << " delta=" << delta << std::endl;

// do a line search
double step_length = 1;
double step_length = 1.0;
do
{
// update of p_unit. The spacedim-th component of transformed
Expand All @@ -630,11 +634,14 @@ namespace internal
f_weighted_norm_square = (df_inverse * f_trial).norm_square();

if (print_iterations_to_deallog)
deallog << " step_length=" << step_length << std::endl
<< " ||f || =" << f.norm() << std::endl
<< " ||f*|| =" << f_trial.norm() << std::endl
<< " ||f*||_A ="
<< std::sqrt(f_weighted_norm_square) << std::endl;
{
deallog << " step_length=" << step_length << std::endl;
if (step_length == 1.0)
deallog << " ||f || =" << f.norm() << std::endl;
deallog << " ||f*|| =" << f_trial.norm() << std::endl
<< " ||f*||_A ="
<< std::sqrt(f_weighted_norm_square) << std::endl;
}

// See if we are making progress with the current step length
// and if not, reduce it by a factor of two and try again.
Expand Down Expand Up @@ -674,7 +681,7 @@ namespace internal
// too small, we give the iteration another try with the
// projection of the initial guess to the unit cell before we give
// up, just like for the negative determinant case.
if (step_length <= 0.05 && try_project_to_unit_cell == false)
if (step_length <= 0.05 && tried_project_to_unit_cell == false)
{
p_unit = GeometryInfo<dim>::project_to_unit_cell(p_unit);
p_real = internal::evaluate_tensor_product_value_and_gradient(
Expand All @@ -686,7 +693,7 @@ namespace internal
f = p_real.first - p;
f_weighted_norm_square = 1.;
last_f_weighted_norm_square = 1;
try_project_to_unit_cell = true;
tried_project_to_unit_cell = true;
continue;
}
else if (step_length <= 0.05)
Expand Down Expand Up @@ -906,10 +913,16 @@ namespace internal
make_array_view(real_support_points));
DerivativeForm<1, spacedim, dim> A_inv =
affine.first.covariant_form().transpose();
coefficients[0] = apply_transformation(A_inv, affine.second);

// The code for evaluation assumes an additional transformation of
// the form (x - normalization_shift) * normalization_length --
// account for this in the definition of the coefficients.
coefficients[0] =
apply_transformation(A_inv, normalization_shift - affine.second);
for (unsigned int d = 0; d < spacedim; ++d)
for (unsigned int e = 0; e < dim; ++e)
coefficients[1 + d][e] = A_inv[e][d];
coefficients[1 + d][e] =
A_inv[e][d] * (1.0 / normalization_length);
is_affine = true;
return;
}
Expand Down Expand Up @@ -960,7 +973,8 @@ namespace internal
Lij_sum += matrix[i][j] * matrix[i][j];
}
AssertThrow(matrix[i][i] - Lij_sum >= 0,
ExcMessage("Matrix not positive definite"));
ExcMessage("Matrix of normal equations not positive "
"definite"));

// Store the inverse in the diagonal since that is the quantity
// needed later in the factorization as well as the forward and
Expand Down Expand Up @@ -1028,10 +1042,30 @@ namespace internal

if (!is_affine)
{
Point<dim, Number> result_affine = result;
for (unsigned int d = 0, c = 0; d < spacedim; ++d)
for (unsigned int e = 0; e <= d; ++e, ++c)
result +=
coefficients[1 + spacedim + c] * (p_scaled[d] * p_scaled[e]);

// Check if the quadratic approximation ends up considerably
// farther outside the unit cell on some or all SIMD lanes than
// the affine approximation - in that case, we switch those
// components back to the affine approximation. Note that the
// quadratic approximation will grow more quickly away from the
// unit cell. We make the selection for each SIMD lane with a
// ternary operation.
const Number distance_to_unit_cell = result.distance_square(
GeometryInfo<dim>::project_to_unit_cell(result));
const Number affine_distance_to_unit_cell =
result_affine.distance_square(
GeometryInfo<dim>::project_to_unit_cell(result_affine));
for (unsigned int d = 0; d < dim; ++d)
result[d] = compare_and_apply_mask<SIMDComparison::greater_than>(
distance_to_unit_cell,
affine_distance_to_unit_cell + 0.5,
result_affine[d],
result[d]);
}
return result;
}
Expand Down
32 changes: 16 additions & 16 deletions tests/mappings/mapping_q_inverse_quadratic_approximation.output
Original file line number Diff line number Diff line change
@@ -1,7 +1,7 @@

DEAL::Testing 2D with point -0.1498310826 0.7854420410
DEAL::Testing on cell 0_1:0 with center 0.8227241336 0.4750000000
DEAL::Exact inverse: 1.679734433 2.008122269
DEAL::Exact inverse: 1.679734433 2.008122260
DEAL::Affine approximation: 1.294871897 6.963559088
DEAL::Inverse quadratic approximation: 1.644896745 2.003800318
DEAL::
Expand All @@ -26,7 +26,7 @@ DEAL::Affine approximation: -0.1371549678 3.966937436
DEAL::Inverse quadratic approximation: -0.2909091648 1.929190316
DEAL::
DEAL::Testing on cell 1_1:1 with center -0.8227241336 -0.4750000000
DEAL::Exact inverse: -1.321431976 2.123538334
DEAL::Exact inverse: -1.321431986 2.123538296
DEAL::Affine approximation: -0.2948718968 13.03644091
DEAL::Inverse quadratic approximation: -0.9655241914 2.275064795
DEAL::
Expand Down Expand Up @@ -68,7 +68,7 @@ DEAL::Affine approximation: -0.1769771533 3.589871026
DEAL::Inverse quadratic approximation: -0.3233858788 1.488970865
DEAL::
DEAL::Testing on cell 1_1:1 with center -0.8227241336 -0.4750000000
DEAL::Exact inverse: -1.321431983 1.631260238
DEAL::Exact inverse: -1.321431986 1.631260249
DEAL::Affine approximation: -0.3445513904 13.22621847
DEAL::Inverse quadratic approximation: -1.067765287 1.852528354
DEAL::
Expand All @@ -78,7 +78,7 @@ DEAL::Affine approximation: -0.2566215243 2.589871026
DEAL::Inverse quadratic approximation: -0.3901601371 0.5348403961
DEAL::
DEAL::Testing on cell 1_1:3 with center -0.7361215932 -0.4250000000
DEAL::Exact inverse: -1.321431983 0.6312603084
DEAL::Exact inverse: -1.321431983 0.6312603085
DEAL::Affine approximation: -0.4439103775 12.22621847
DEAL::Inverse quadratic approximation: -1.324499457 0.9907864881
DEAL::
Expand Down Expand Up @@ -115,7 +115,7 @@ DEAL::Affine approximation: -0.3942308840 13.41599603
DEAL::Inverse quadratic approximation: -1.171258833 1.403869886
DEAL::
DEAL::Testing on cell 1_1:2 with center -0.7361215932 0.4250000000
DEAL::Exact inverse: -0.3199797942 0.004618771511
DEAL::Exact inverse: -0.3199797942 0.004618771512
DEAL::Affine approximation: -0.3011286728 2.212804616
DEAL::Inverse quadratic approximation: -0.4188102773 0.01487792182
DEAL::
Expand All @@ -137,7 +137,7 @@ DEAL::Affine approximation: 0.6872888532 -0.7700353780
DEAL::Inverse quadratic approximation: 0.6735562329 0.5114552797
DEAL::
DEAL::Testing on cell 0_1:2 with center 0.7361215932 0.4250000000
DEAL::Exact inverse: 1.679734434 -0.4903548181
DEAL::Exact inverse: 1.679734434 -0.4903548182
DEAL::Affine approximation: 1.554958657 5.394226417
DEAL::Inverse quadratic approximation: 1.973382985 -0.4493850397
DEAL::
Expand All @@ -162,14 +162,14 @@ DEAL::Affine approximation: -0.3456358213 1.835738205
DEAL::Inverse quadratic approximation: -0.4448257264 -0.5335632508
DEAL::
DEAL::Testing on cell 1_1:3 with center -0.7361215932 -0.4250000000
DEAL::Exact inverse: -1.321431982 -0.3532973076
DEAL::Exact inverse: -1.321431983 -0.3532973076
DEAL::Affine approximation: -0.5549586572 12.60577358
DEAL::Inverse quadratic approximation: -1.570577755 -0.03655836469
DEAL::Inverse quadratic approximation: -1.570577755 -0.03655836468
DEAL::
DEAL::
DEAL::Testing 2D with point -0.1872888532 0.9818025513
DEAL::Testing on cell 0_1:0 with center 0.8227241336 0.4750000000
DEAL::Exact inverse: 1.679734435 0.01015282302
DEAL::Exact inverse: 1.679734435 0.01015282301
DEAL::Affine approximation: 1.493589871 6.204448860
DEAL::Inverse quadratic approximation: 1.881022893 0.08910626957
DEAL::
Expand All @@ -179,14 +179,14 @@ DEAL::Affine approximation: 0.6971461613 -1.336879345
DEAL::Inverse quadratic approximation: 0.6752680942 0.001010067356
DEAL::
DEAL::Testing on cell 0_1:2 with center 0.7361215932 0.4250000000
DEAL::Exact inverse: 1.679734433 -0.9898471766
DEAL::Exact inverse: 1.679734433 -0.9898471760
DEAL::Affine approximation: 1.610482797 5.204448860
DEAL::Inverse quadratic approximation: 2.034386667 -1.028496543
DEAL::
DEAL::Testing on cell 0_1:3 with center 5.204748896e-17 0.8500000000
DEAL::Exact inverse: 0.6800000317 -0.9950656588
DEAL::Affine approximation: 0.7203398273 -2.336879345
DEAL::Inverse quadratic approximation: 0.6740153620 -1.104288093
DEAL::Inverse quadratic approximation: 0.8115480065 1.507109783
DEAL::
DEAL::Testing on cell 1_1:0 with center -0.8227241336 0.4750000000
DEAL::Exact inverse: -0.3199797942 0.005131968346
Expand All @@ -199,12 +199,12 @@ DEAL::Affine approximation: -0.4935898711 13.79555114
DEAL::Inverse quadratic approximation: -1.382003278 0.4281868663
DEAL::
DEAL::Testing on cell 1_1:2 with center -0.7361215932 0.4250000000
DEAL::Exact inverse: -0.3199797940 -0.9948680322
DEAL::Exact inverse: -0.3199797941 -0.9948680324
DEAL::Affine approximation: -0.3901429698 1.458671795
DEAL::Inverse quadratic approximation: -0.4682064843 -1.110483122
DEAL::Inverse quadratic approximation: -0.4269023613 -0.4705876412
DEAL::
DEAL::Testing on cell 1_1:3 with center -0.7361215932 -0.4250000000
DEAL::Exact inverse: -1.321431978 -0.8455762513
DEAL::Exact inverse: -1.321431974 -0.8455762606
DEAL::Affine approximation: -0.6104827971 12.79555114
DEAL::Inverse quadratic approximation: -1.696098323 -0.5939421483
DEAL::
Expand Down Expand Up @@ -295,7 +295,7 @@ DEAL::
DEAL::
DEAL::Testing 3D with point 0.05085652866 0.06999800658 0.8623369432
DEAL::Testing on cell 0_0: with center -9.723960037e-18 -9.723960037e-18 -0.9000000000
DEAL::Exact inverse: 0.4633821073 0.4496318287 9.333474387
DEAL::Exact inverse: 0.4633821073 0.4496318287 9.333474388
DEAL::Affine approximation: 0.5489367175 0.5673556132 12.46805699
DEAL::Inverse quadratic approximation: 0.6217110739 0.6675209216 0.8911022773
DEAL::
Expand Down Expand Up @@ -331,7 +331,7 @@ DEAL::
DEAL::
DEAL::Testing 3D with point 0.05672458966 0.07807469965 0.9618373598
DEAL::Testing on cell 0_0: with center -9.723960037e-18 -9.723960037e-18 -0.9000000000
DEAL::Exact inverse: 0.4633821073 0.4496318287 9.833490663
DEAL::Exact inverse: 0.4633821073 0.4496318287 9.833490662
DEAL::Affine approximation: 0.5545832618 0.5751274148 13.32975588
DEAL::Inverse quadratic approximation: 0.6416309076 0.6949382206 0.4108019385
DEAL::
Expand Down
75 changes: 75 additions & 0 deletions tests/mappings/mapping_q_points_real_to_unit_02.cc
Original file line number Diff line number Diff line change
@@ -0,0 +1,75 @@
// ---------------------------------------------------------------------
//
// Copyright (C) 2020 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.md at
// the top level directory of deal.II.
//
// ---------------------------------------------------------------------


// check MappingQ::transform_real_to_unit_points on MappingQ1 with a point
// that would at some point return an invalid answer due to a wrongly computed
// affine approximation, using a point extracted from
// tests/particles/particle_handler_21.

#include <deal.II/fe/mapping_q.h>

#include <deal.II/grid/grid_generator.h>
#include <deal.II/grid/tria.h>

#include "../tests.h"


template <int dim>
void
test()
{
Triangulation<dim> tria;
GridGenerator::hyper_ball(tria, Point<dim>(), 1.0);
tria.refine_global(3);

// pick out a particular cell
const typename Triangulation<dim>::cell_iterator coarse_cell(&tria, 0, 5);
const auto cell = coarse_cell->child(0)->child(5)->child(7);

std::vector<Point<dim>> real_points{
Point<3>(-0.0489699, -0.859053, -0.0489699)};
std::vector<Point<dim>> unit_points(real_points.size());

for (unsigned int degree = 1; degree < 4; ++degree)
{
MappingQ<dim> mapping(degree);
mapping.transform_points_real_to_unit_cell(cell,
real_points,
unit_points);

deallog << "Transform on cell with center: " << cell->center(true)
<< " with mapping degree " << degree << std::endl;
deallog << "Combined transform " << real_points[0] << " gives "
<< unit_points[0] << std::endl;

deallog << "Classic transform " << real_points[0] << " gives "
<< mapping.transform_real_to_unit_cell(cell, real_points[0])
<< std::endl;
}

deallog << "OK" << std::endl;
}



int
main()
{
initlog();
deallog << std::setprecision(10);

test<3>();
}
11 changes: 11 additions & 0 deletions tests/mappings/mapping_q_points_real_to_unit_02.output
Original file line number Diff line number Diff line change
@@ -0,0 +1,11 @@

DEAL::Transform on cell with center: -0.08086235035 -0.8200448543 -0.08086235035 with mapping degree 1
DEAL::Combined transform -0.04896990000 -0.8590530000 -0.04896990000 gives 0.7082205181 0.2510034894 0.7082205181
DEAL::Classic transform -0.04896990000 -0.8590530000 -0.04896990000 gives 0.7082205181 0.2510034894 0.7082205181
DEAL::Transform on cell with center: -0.08086235035 -0.8200448543 -0.08086235035 with mapping degree 2
DEAL::Combined transform -0.04896990000 -0.8590530000 -0.04896990000 gives 0.7082205181 0.2510034894 0.7082205181
DEAL::Classic transform -0.04896990000 -0.8590530000 -0.04896990000 gives 0.7082205181 0.2510034894 0.7082205181
DEAL::Transform on cell with center: -0.08086235035 -0.8200448543 -0.08086235035 with mapping degree 3
DEAL::Combined transform -0.04896990000 -0.8590530000 -0.04896990000 gives 0.7082205181 0.2510034894 0.7082205181
DEAL::Classic transform -0.04896990000 -0.8590530000 -0.04896990000 gives 0.7082205181 0.2510034894 0.7082205181
DEAL::OK

0 comments on commit 4e8af90

Please sign in to comment.