]> https://gitweb.dealii.org/ - dealii.git/commitdiff
MappingFE::transform_real_to_unit_cell()
authorLuca Heltai <luca.heltai@sissa.it>
Wed, 13 Jan 2021 17:28:30 +0000 (18:28 +0100)
committerLuca Heltai <luca.heltai@sissa.it>
Thu, 14 Jan 2021 15:59:41 +0000 (16:59 +0100)
doc/news/changes/minor/20210113LucaHeltai [new file with mode: 0644]
include/deal.II/fe/mapping_q_internal.h
source/fe/mapping_fe.cc
tests/simplex/mapping_transformations_01.cc [new file with mode: 0644]
tests/simplex/mapping_transformations_01.with_simplex_support=on.output [new file with mode: 0644]

diff --git a/doc/news/changes/minor/20210113LucaHeltai b/doc/news/changes/minor/20210113LucaHeltai
new file mode 100644 (file)
index 0000000..dbd3d45
--- /dev/null
@@ -0,0 +1,3 @@
+New: Implemented MappingFE::transform_real_to_unit_cell().
+<br>
+(Luca Heltai, 2021/01/13)
index d4d86f90802ad25a4b476711e08ff658f68ecf50..e747f747a61e37a304236516e78f8b9a4ef4ecdf 100644 (file)
@@ -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
index 3f7148dbca07893aad0e98d37c88592c8311fd86..7d605bde60a1f7d97e87742ad9c6ba188972f801 100644 (file)
@@ -909,12 +909,84 @@ MappingFE<dim, spacedim>::transform_real_to_unit_cell(
   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;
 }
 
 
@@ -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<dim, spacedim>::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<dim, spacedim>::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 (file)
index 0000000..830fa8a
--- /dev/null
@@ -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 <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;
+      }
+}
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 (file)
index 0000000..eef90cd
--- /dev/null
@@ -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

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.