const typename Triangulation<dim, spacedim>::cell_iterator &cell,
const Point<spacedim> & p) const
{
- Assert(false, StandardExceptions::ExcNotImplemented());
+ const std::vector<Point<spacedim>> support_points =
+ this->compute_mapping_support_points(cell);
+
+ const double eps = 1.e-12 * cell->diameter();
+ const unsigned int loop_limit = 10;
+
+ Point<dim> 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<spacedim> 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<dim>(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<dim, spacedim>::ExcTransformationFailed()));
- return {};
+ return p_unit;
}
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.
*/
// 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 <
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)
// 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;
}
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
// 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()) ||
// 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()) ||
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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 <deal.II/base/patterns.h>
+
+#include <deal.II/fe/mapping_fe.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/tria_accessor.h>
+#include <deal.II/grid/tria_iterator.h>
+
+#include <deal.II/simplex/fe_lib.h>
+
+#include <fstream>
+
+#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<Point<2>> 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;
+ }
+}
--- /dev/null
+
+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
projects / dealii.git / commitdiff