From c36043106a8235e72f32e052c849162b1f7db3c6 Mon Sep 17 00:00:00 2001 From: Luca Heltai Date: Wed, 13 Jan 2021 18:28:30 +0100 Subject: [PATCH] MappingFE::transform_real_to_unit_cell() --- doc/news/changes/minor/20210113LucaHeltai | 3 + include/deal.II/fe/mapping_q_internal.h | 2 +- source/fe/mapping_fe.cc | 116 ++++++++++++++---- tests/simplex/mapping_transformations_01.cc | 86 +++++++++++++ ...rmations_01.with_simplex_support=on.output | 43 +++++++ 5 files changed, 228 insertions(+), 22 deletions(-) create mode 100644 doc/news/changes/minor/20210113LucaHeltai create mode 100644 tests/simplex/mapping_transformations_01.cc create mode 100644 tests/simplex/mapping_transformations_01.with_simplex_support=on.output diff --git a/doc/news/changes/minor/20210113LucaHeltai b/doc/news/changes/minor/20210113LucaHeltai new file mode 100644 index 0000000000..dbd3d4510c --- /dev/null +++ b/doc/news/changes/minor/20210113LucaHeltai @@ -0,0 +1,3 @@ +New: Implemented MappingFE::transform_real_to_unit_cell(). +
+(Luca Heltai, 2021/01/13) diff --git a/include/deal.II/fe/mapping_q_internal.h b/include/deal.II/fe/mapping_q_internal.h index d4d86f9080..e747f747a6 100644 --- a/include/deal.II/fe/mapping_q_internal.h +++ b/include/deal.II/fe/mapping_q_internal.h @@ -531,7 +531,7 @@ namespace internal // p(x) - p = p(x) - p(x*) // = -grad p(x) * (x*-x) + higher order terms // This suggest to measure with a norm that corresponds to - // A = {[grad p(x]^T [grad p(x)]}^{-1} + // A = {[grad p(x)]^T [grad p(x)]}^{-1} // because then // \| p(x) - p \|_A \approx \| x - x* \| // Consequently, we will try to enforce that diff --git a/source/fe/mapping_fe.cc b/source/fe/mapping_fe.cc index 3f7148dbca..7d605bde60 100644 --- a/source/fe/mapping_fe.cc +++ b/source/fe/mapping_fe.cc @@ -909,12 +909,84 @@ MappingFE::transform_real_to_unit_cell( const typename Triangulation::cell_iterator &cell, const Point & p) const { - Assert(false, StandardExceptions::ExcNotImplemented()); + const std::vector> support_points = + this->compute_mapping_support_points(cell); + + const double eps = 1.e-12 * cell->diameter(); + const unsigned int loop_limit = 10; + + Point p_unit; + + Tensor<1, dim> grad_FT_residual; + + unsigned int loop = 0; + + // This loop solves the following problem: + // grad_F^T residual + (grad_F^T grad_F + grad_F^T hesss_F dp) dp = 0 + // where the term + // (grad_F^T hess_F dp) is approximated by (hess_F * residual) + // This is basically a second order approximation of Newton method, where the + // Jacobian is corrected with a higher order term coming from the hessian. + do + { + Point mapped_point; + + // Transpose of the gradient map + DerivativeForm<1, spacedim, dim> grad_FT; + Tensor<2, dim> corrected_metric_tensor; + + // [TODO] + // When 2nd derivatives are implemented, we'll uncomment this + // DerivativeForm<2, spacedim, dim> hess_FT; + for (unsigned int i = 0; i < this->fe->n_dofs_per_cell(); ++i) + { + mapped_point += support_points[i] * this->fe->shape_value(i, p_unit); + const auto grad_F_i = this->fe->shape_grad(i, p_unit); + // [TODO] + // When 2nd derivatives are implemented, we'll uncomment this + // auto hessian_k = this->fe->shape_grad_grad(i, p_unit); + for (unsigned int j = 0; j < dim; ++j) + { + grad_FT[j] += grad_F_i[j] * support_points[i]; + // [TODO] + // When 2nd derivatives are implemented, we'll uncomment this + // hess_FT[j] += hessian_k[j] * support_points[i]; + } + } + + // Residual + auto residual = p - mapped_point; + // Project the residual on the reference coordinate system + // to compute the error, and to filter components orthogonal to the + // manifold, and compute a 2nd order correction of the metric tensor + grad_FT_residual = apply_transformation(grad_FT, residual); + + // Do not invert nor compute the metric if not necessary. + if (grad_FT_residual.norm() <= eps) + break; + + // Now compute the (corrected) metric tensor + corrected_metric_tensor = -apply_transformation(grad_FT, grad_FT); + + // [TODO] + // When 2nd derivatives are implemented, we'll uncomment this + // corrected_metric_tensor += apply_transformation(hess_FT, residual); + + auto g_inverse = invert(corrected_metric_tensor); + p_unit -= Point(g_inverse * grad_FT_residual); + + ++loop; + } + while (loop < loop_limit); - (void)cell; - (void)p; + // Here we check that in the last execution of while the first + // condition was already wrong, meaning the residual was below + // eps. Only if the first condition failed, loop will have been + // increased and tested, and thus have reached the limit. + AssertThrow(loop < loop_limit, + (typename Mapping::ExcTransformationFailed())); - return {}; + return p_unit; } @@ -1250,11 +1322,12 @@ namespace internal namespace { /** - * Depending on what information is called for in the update flags of the + * Depending on what information is called for in the update flags of + * the * @p data object, compute the various pieces of information that is * required by the fill_fe_face_values() and fill_fe_subface_values() - * functions. This function simply unifies the work that would be done by - * those two functions. + * functions. This function simply unifies the work that would be done + * by those two functions. * * The resulting data is put into the @p output_data argument. */ @@ -1292,9 +1365,9 @@ namespace internal // first compute some common data that is used for evaluating // all of the flags below - // map the unit tangentials to the real cell. checking for d!=dim-1 - // eliminates compiler warnings regarding unsigned int expressions < - // 0. + // map the unit tangentials to the real cell. checking for + // d!=dim-1 eliminates compiler warnings regarding unsigned int + // expressions < 0. for (unsigned int d = 0; d != dim - 1; ++d) { Assert(face_no + cell->n_faces() * d < @@ -1315,8 +1388,9 @@ namespace internal if (update_flags & update_boundary_forms) { - // if dim==spacedim, we can use the unit tangentials to compute - // the boundary form by simply taking the cross product + // if dim==spacedim, we can use the unit tangentials to + // compute the boundary form by simply taking the cross + // product if (dim == spacedim) { for (unsigned int i = 0; i < n_q_points; ++i) @@ -1326,8 +1400,8 @@ namespace internal // in 1d, we don't have access to any of the // data.aux fields (because it has only dim-1 // components), but we can still compute the - // boundary form by simply looking at the number of - // the face + // boundary form by simply looking at the number + // of the face output_data.boundary_forms[i][0] = (face_no == 0 ? -1 : +1); break; @@ -1345,9 +1419,9 @@ namespace internal } else //(dim < spacedim) { - // in the codim-one case, the boundary form results from the - // cross product of all the face tangential vectors and the - // cell normal vector + // in the codim-one case, the boundary form results from + // the cross product of all the face tangential vectors + // and the cell normal vector // // to compute the cell normal, use the same method used in // fill_fe_values for cells above @@ -1519,8 +1593,8 @@ MappingFE::fill_fe_face_values( // if necessary, recompute the support points of the transformation of this // cell (note that we need to first check the triangulation pointer, since - // otherwise the second test might trigger an exception if the triangulations - // are not the same) + // otherwise the second test might trigger an exception if the + // triangulations are not the same) if ((data.mapping_support_points.size() == 0) || (&cell->get_triangulation() != &data.cell_of_current_support_points->get_triangulation()) || @@ -1566,8 +1640,8 @@ MappingFE::fill_fe_subface_values( // if necessary, recompute the support points of the transformation of this // cell (note that we need to first check the triangulation pointer, since - // otherwise the second test might trigger an exception if the triangulations - // are not the same) + // otherwise the second test might trigger an exception if the + // triangulations are not the same) if ((data.mapping_support_points.size() == 0) || (&cell->get_triangulation() != &data.cell_of_current_support_points->get_triangulation()) || diff --git a/tests/simplex/mapping_transformations_01.cc b/tests/simplex/mapping_transformations_01.cc new file mode 100644 index 0000000000..830fa8a4e1 --- /dev/null +++ b/tests/simplex/mapping_transformations_01.cc @@ -0,0 +1,86 @@ +// --------------------------------------------------------------------- +// +// 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. +// +// --------------------------------------------------------------------- + + +// Create a simplex mesh, and transform back and forth a point in each cell + +#include + +#include + +#include +#include +#include +#include + +#include + +#include + +#include "../tests.h" + +using namespace dealii; + +void make_grid(Triangulation<2> &triangulation) +{ + Triangulation<2> triangulation_temp; + + const Point<2> center(1, 0); + const double inner_radius = 0.5, outer_radius = 1.0; + GridGenerator::hyper_shell( + triangulation_temp, center, inner_radius, outer_radius, 5); + + GridGenerator::convert_hypercube_to_simplex_mesh(triangulation_temp, + triangulation); +} + +int +main() +{ + initlog(); + + Triangulation<2> triangulation; + make_grid(triangulation); + + MappingFE<2> mapping(Simplex::FE_P<2>(1)); + + unsigned int n_points = 1; + + std::vector> points; + do + { + auto p = random_point<2>(); + auto sum = p[0] + p[1]; + if (sum <= 1) + points.push_back(p); + } + while (points.size() < n_points); + + deallog << "Transforming " << n_points << " from reference to real." + << std::endl + << "Points: " << Patterns::Tools::to_string(points) << std::endl; + + for (const auto &cell : triangulation.active_cell_iterators()) + for (const auto &p : points) + { + const auto real_p = mapping.transform_unit_to_real_cell(cell, p); + const auto pull_back_p = + mapping.transform_real_to_unit_cell(cell, real_p); + deallog << "F(" << p << ")=" << real_p << " --- " + << "F^-1(" << real_p << ")=" << pull_back_p << std::endl; + if (p.distance(pull_back_p) > 1e-10) + deallog << "ERROR: " << p.distance(pull_back_p) << std::endl; + } +} diff --git a/tests/simplex/mapping_transformations_01.with_simplex_support=on.output b/tests/simplex/mapping_transformations_01.with_simplex_support=on.output new file mode 100644 index 0000000000..eef90cd837 --- /dev/null +++ b/tests/simplex/mapping_transformations_01.with_simplex_support=on.output @@ -0,0 +1,43 @@ + +DEAL::Transforming 1 from reference to real. +DEAL::Points: 0.277775, 0.55397 +DEAL::F(0.277775 0.553970)=1.80846 0.163272 --- F^-1(1.80846 0.163272)=0.277775 0.553970 +DEAL::F(0.277775 0.553970)=1.75859 0.399789 --- F^-1(1.75859 0.399789)=0.277775 0.553970 +DEAL::F(0.277775 0.553970)=1.45521 0.632588 --- F^-1(1.45521 0.632588)=0.277775 0.553970 +DEAL::F(0.277775 0.553970)=1.56454 0.683770 --- F^-1(1.56454 0.683770)=0.277775 0.553970 +DEAL::F(0.277775 0.553970)=1.51654 0.162808 --- F^-1(1.51654 0.162808)=0.277775 0.553970 +DEAL::F(0.277775 0.553970)=1.62993 0.155809 --- F^-1(1.62993 0.155809)=0.277775 0.553970 +DEAL::F(0.277775 0.553970)=1.39055 0.435116 --- F^-1(1.39055 0.435116)=0.277775 0.553970 +DEAL::F(0.277775 0.553970)=1.26675 0.556789 --- F^-1(1.26675 0.556789)=0.277775 0.553970 +DEAL::F(0.277775 0.553970)=1.09455 0.819342 --- F^-1(1.09455 0.819342)=0.277775 0.553970 +DEAL::F(0.277775 0.553970)=0.854196 0.845007 --- F^-1(0.854196 0.845007)=0.277775 0.553970 +DEAL::F(0.277775 0.553970)=0.539039 0.628406 --- F^-1(0.539039 0.628406)=0.277775 0.553970 +DEAL::F(0.277775 0.553970)=0.524149 0.748209 --- F^-1(0.524149 0.748209)=0.277775 0.553970 +DEAL::F(0.277775 0.553970)=1.00478 0.541573 --- F^-1(1.00478 0.541573)=0.277775 0.553970 +DEAL::F(0.277775 0.553970)=1.04648 0.647248 --- F^-1(1.04648 0.647248)=0.277775 0.553970 +DEAL::F(0.277775 0.553970)=0.706868 0.505898 --- F^-1(0.706868 0.505898)=0.277775 0.553970 +DEAL::F(0.277775 0.553970)=0.552892 0.425750 --- F^-1(0.552892 0.425750)=0.277775 0.553970 +DEAL::F(0.277775 0.553970)=0.249976 0.343110 --- F^-1(0.249976 0.343110)=0.277775 0.553970 +DEAL::F(0.277775 0.553970)=0.151295 0.122454 --- F^-1(0.151295 0.122454)=0.277775 0.553970 +DEAL::F(0.277775 0.553970)=0.259905 -0.244212 --- F^-1(0.259905 -0.244212)=0.277775 0.553970 +DEAL::F(0.277775 0.553970)=0.141365 -0.221352 --- F^-1(0.141365 -0.221352)=0.277775 0.553970 +DEAL::F(0.277775 0.553970)=0.486411 0.171903 --- F^-1(0.486411 0.171903)=0.277775 0.553970 +DEAL::F(0.277775 0.553970)=0.398793 0.244212 --- F^-1(0.398793 0.244212)=0.277775 0.553970 +DEAL::F(0.277775 0.553970)=0.428280 -0.122454 --- F^-1(0.428280 -0.122454)=0.277775 0.553970 +DEAL::F(0.277775 0.553970)=0.456924 -0.293661 --- F^-1(0.456924 -0.293661)=0.277775 0.553970 +DEAL::F(0.277775 0.553970)=0.441913 -0.607289 --- F^-1(0.441913 -0.607289)=0.277775 0.553970 +DEAL::F(0.277775 0.553970)=0.621275 -0.769326 --- F^-1(0.621275 -0.769326)=0.277775 0.553970 +DEAL::F(0.277775 0.553970)=1.00356 -0.779337 --- F^-1(1.00356 -0.779337)=0.277775 0.553970 +DEAL::F(0.277775 0.553970)=0.945185 -0.885012 --- F^-1(0.945185 -0.885012)=0.277775 0.553970 +DEAL::F(0.277775 0.553970)=0.677803 -0.435331 --- F^-1(0.677803 -0.435331)=0.277775 0.553970 +DEAL::F(0.277775 0.553970)=0.581958 -0.496317 --- F^-1(0.581958 -0.496317)=0.277775 0.553970 +DEAL::F(0.277775 0.553970)=0.939789 -0.581578 --- F^-1(0.939789 -0.581578)=0.277775 0.553970 +DEAL::F(0.277775 0.553970)=1.11147 -0.607242 --- F^-1(1.11147 -0.607242)=0.277775 0.553970 +DEAL::F(0.277775 0.553970)=1.40511 -0.718435 --- F^-1(1.40511 -0.718435)=0.277775 0.553970 +DEAL::F(0.277775 0.553970)=1.61464 -0.597923 --- F^-1(1.61464 -0.597923)=0.277775 0.553970 +DEAL::F(0.277775 0.553970)=1.74229 -0.237445 --- F^-1(1.74229 -0.237445)=0.277775 0.553970 +DEAL::F(0.277775 0.553970)=1.82476 -0.325615 --- F^-1(1.82476 -0.325615)=0.277775 0.553970 +DEAL::F(0.277775 0.553970)=1.31446 -0.440952 --- F^-1(1.31446 -0.440952)=0.277775 0.553970 +DEAL::F(0.277775 0.553970)=1.34284 -0.550952 --- F^-1(1.34284 -0.550952)=0.277775 0.553970 +DEAL::F(0.277775 0.553970)=1.53451 -0.236981 --- F^-1(1.53451 -0.236981)=0.277775 0.553970 +DEAL::F(0.277775 0.553970)=1.61197 -0.0816359 --- F^-1(1.61197 -0.0816359)=0.277775 0.553970 -- 2.39.5