]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Move the implementation of transform_r_to_u_cell from MappingQ1 to MappingQGeneric.
authorWolfgang Bangerth <bangerth@math.tamu.edu>
Sat, 19 Sep 2015 13:03:15 +0000 (08:03 -0500)
committerWolfgang Bangerth <bangerth@math.tamu.edu>
Thu, 1 Oct 2015 20:16:34 +0000 (15:16 -0500)
include/deal.II/fe/mapping_q1.h
include/deal.II/fe/mapping_q_generic.h
source/fe/mapping_q1.cc
source/fe/mapping_q_generic.cc

index ba05c5cfc7dba410198f86fd02fcd8cd64450b80..33371651a276824c1d15ddce24ab3c1bb0d63c60 100644 (file)
@@ -63,21 +63,6 @@ public:
   // for documentation, see the Mapping base class
   virtual
   MappingQ1<dim,spacedim> *clone () const;
-
-  /**
-   * @name Mapping points between reference and real cells
-   * @{
-   */
-
-  // for documentation, see the Mapping base class
-  virtual
-  Point<dim>
-  transform_real_to_unit_cell (const typename Triangulation<dim,spacedim>::cell_iterator &cell,
-                               const Point<spacedim>                            &p) const;
-
-  /**
-   * @}
-   */
 };
 
 
index 8b876feebff73411029f5ae825cca379bc84c2f9..d3fa8a4c21e7af1a748c893e9f8a2339822bd143 100644 (file)
@@ -96,6 +96,12 @@ public:
   transform_unit_to_real_cell (const typename Triangulation<dim,spacedim>::cell_iterator &cell,
                                const Point<dim>                                 &p) const;
 
+  // for documentation, see the Mapping base class
+  virtual
+  Point<dim>
+  transform_real_to_unit_cell (const typename Triangulation<dim,spacedim>::cell_iterator &cell,
+                               const Point<spacedim>                            &p) const;
+
   /**
    * @}
    */
index 54001d53b2af58b7ac3019db56d990c8a454ddf5..e61f37e3399723d9035fad20260a1ff3f2d0a08d 100644 (file)
@@ -48,401 +48,6 @@ MappingQ1<dim,spacedim>::MappingQ1 ()
 
 
 
-namespace internal
-{
-  namespace MappingQ1
-  {
-    // These are left as templates on the spatial dimension (even though dim
-    // == spacedim must be true for them to make sense) because templates are
-    // expanded before the compiler eliminates code due to the 'if (dim ==
-    // spacedim)' statement (see the body of the general
-    // transform_real_to_unit_cell).
-    template<int spacedim>
-    Point<1>
-    transform_real_to_unit_cell
-    (const std_cxx11::array<Point<spacedim>, GeometryInfo<1>::vertices_per_cell> &vertices,
-     const Point<spacedim> &p)
-    {
-      Assert(spacedim == 1, ExcInternalError());
-      return Point<1>((p[0] - vertices[0](0))/(vertices[1](0) - vertices[0](0)));
-    }
-
-
-
-    template<int spacedim>
-    Point<2>
-    transform_real_to_unit_cell
-    (const std_cxx11::array<Point<spacedim>, GeometryInfo<2>::vertices_per_cell> &vertices,
-     const Point<spacedim> &p)
-    {
-      Assert(spacedim == 2, ExcInternalError());
-      const double x = p(0);
-      const double y = p(1);
-
-      const double x0 = vertices[0](0);
-      const double x1 = vertices[1](0);
-      const double x2 = vertices[2](0);
-      const double x3 = vertices[3](0);
-
-      const double y0 = vertices[0](1);
-      const double y1 = vertices[1](1);
-      const double y2 = vertices[2](1);
-      const double y3 = vertices[3](1);
-
-      const double a = (x1 - x3)*(y0 - y2) - (x0 - x2)*(y1 - y3);
-      const double b = -(x0 - x1 - x2 + x3)*y + (x - 2*x1 + x3)*y0 - (x - 2*x0 + x2)*y1
-                       - (x - x1)*y2 + (x - x0)*y3;
-      const double c = (x0 - x1)*y - (x - x1)*y0 + (x - x0)*y1;
-
-      const double discriminant = b*b - 4*a*c;
-      // exit if the point is not in the cell (this is the only case where the
-      // discriminant is negative)
-      if (discriminant < 0.0)
-        {
-          AssertThrow (false,
-                       (typename Mapping<spacedim,spacedim>::ExcTransformationFailed()));
-        }
-
-      double eta1;
-      double eta2;
-      // special case #1: if a is zero, then use the linear formula
-      if (a == 0.0 && b != 0.0)
-        {
-          eta1 = -c/b;
-          eta2 = -c/b;
-        }
-      // special case #2: if c is very small:
-      else if (std::abs(c/b) < 1e-12)
-        {
-          eta1 = (-b - std::sqrt(discriminant)) / (2*a);
-          eta2 = (-b + std::sqrt(discriminant)) / (2*a);
-        }
-      // finally, use the numerically stable version of the quadratic formula:
-      else
-        {
-          eta1 = 2*c / (-b - std::sqrt(discriminant));
-          eta2 = 2*c / (-b + std::sqrt(discriminant));
-        }
-      // pick the one closer to the center of the cell.
-      const double eta = (std::abs(eta1 - 0.5) < std::abs(eta2 - 0.5)) ? eta1 : eta2;
-
-      /*
-       * There are two ways to compute xi from eta, but either one may have a
-       * zero denominator.
-       */
-      const double subexpr0 = -eta*x2 + x0*(eta - 1);
-      const double xi_denominator0 = eta*x3 - x1*(eta - 1) + subexpr0;
-      const double max_x = std::max(std::max(std::abs(x0), std::abs(x1)),
-                                    std::max(std::abs(x2), std::abs(x3)));
-
-      if (std::abs(xi_denominator0) > 1e-10*max_x)
-        {
-          const double xi = (x + subexpr0)/xi_denominator0;
-          return Point<2>(xi, eta);
-        }
-      else
-        {
-          const double max_y = std::max(std::max(std::abs(y0), std::abs(y1)),
-                                        std::max(std::abs(y2), std::abs(y3)));
-          const double subexpr1 = -eta*y2 + y0*(eta - 1);
-          const double xi_denominator1 = eta*y3 - y1*(eta - 1) + subexpr1;
-          if (std::abs(xi_denominator1) > 1e-10*max_y)
-            {
-              const double xi = (subexpr1 + y)/xi_denominator1;
-              return Point<2>(xi, eta);
-            }
-          else // give up and try Newton iteration
-            {
-              AssertThrow (false,
-                           (typename Mapping<spacedim,spacedim>::ExcTransformationFailed()));
-            }
-        }
-      // bogus return to placate compiler. It should not be possible to get
-      // here.
-      Assert(false, ExcInternalError());
-      return Point<2>(std::numeric_limits<double>::quiet_NaN(),
-                      std::numeric_limits<double>::quiet_NaN());
-    }
-
-
-
-    template<int spacedim>
-    Point<3>
-    transform_real_to_unit_cell
-    (const std_cxx11::array<Point<spacedim>, GeometryInfo<3>::vertices_per_cell> &/*vertices*/,
-     const Point<spacedim> &/*p*/)
-    {
-      // It should not be possible to get here
-      Assert(false, ExcInternalError());
-      return Point<3>();
-    }
-  }
-}
-
-
-
-/**
- * Compute an initial guess to pass to the Newton method in
- * transform_real_to_unit_cell.  For the initial guess we proceed in the
- * following way:
- * <ul>
- * <li> find the least square dim-dimensional plane approximating the cell
- * vertices, i.e. we find an affine map A x_hat + b from the reference cell
- * to the real space.
- * <li> Solve the equation A x_hat + b = p for x_hat
- * <li> This x_hat is the initial solution used for the Newton Method.
- * </ul>
- *
- * @note if dim<spacedim we first project p onto the plane.
- *
- * @note if dim==1 (for any spacedim) the initial guess is the exact
- * solution and no Newton iteration is needed.
- *
- * Some details about how we compute the least square plane. We look
- * for a spacedim x (dim + 1) matrix X such that X * M = Y where M is
- * a (dim+1) x n_vertices matrix and Y a spacedim x n_vertices.  And:
- * The i-th column of M is unit_vertex[i] and the last row all
- * 1's. The i-th column of Y is real_vertex[i].  If we split X=[A|b],
- * the least square approx is A x_hat+b Classically X = Y * (M^t (M
- * M^t)^{-1}) Let K = M^t * (M M^t)^{-1} = [KA Kb] this can be
- * precomputed, and that is exactly what we do.  Finally A = Y*KA and
- * b = Y*Kb.
- */
-namespace
-{
-  template <int dim>
-  struct TransformR2UInitialGuess
-  {
-    static const double KA[GeometryInfo<dim>::vertices_per_cell][dim];
-    static const double Kb[GeometryInfo<dim>::vertices_per_cell];
-  };
-
-
-  /*
-    Octave code:
-    M=[0 1; 1 1];
-    K1 = transpose(M) * inverse (M*transpose(M));
-    printf ("{%f, %f},\n", K1' );
-  */
-  template <>
-  const double
-  TransformR2UInitialGuess<1>::
-  KA[GeometryInfo<1>::vertices_per_cell][1] =
-  {
-    {-1.000000},
-    {1.000000}
-  };
-
-  template <>
-  const double
-  TransformR2UInitialGuess<1>::
-  Kb[GeometryInfo<1>::vertices_per_cell] = {1.000000, 0.000000};
-
-
-  /*
-    Octave code:
-    M=[0 1 0 1;0 0 1 1;1 1 1 1];
-    K2 = transpose(M) * inverse (M*transpose(M));
-    printf ("{%f, %f, %f},\n", K2' );
-  */
-  template <>
-  const double
-  TransformR2UInitialGuess<2>::
-  KA[GeometryInfo<2>::vertices_per_cell][2] =
-  {
-    {-0.500000, -0.500000},
-    { 0.500000, -0.500000},
-    {-0.500000,  0.500000},
-    { 0.500000,  0.500000}
-  };
-
-  /*
-    Octave code:
-    M=[0 1 0 1 0 1 0 1;0 0 1 1 0 0 1 1; 0 0 0 0 1 1 1 1; 1 1 1 1 1 1 1 1];
-    K3 = transpose(M) * inverse (M*transpose(M))
-    printf ("{%f, %f, %f, %f},\n", K3' );
-  */
-  template <>
-  const double
-  TransformR2UInitialGuess<2>::
-  Kb[GeometryInfo<2>::vertices_per_cell] =
-  {0.750000,0.250000,0.250000,-0.250000 };
-
-
-  template <>
-  const double
-  TransformR2UInitialGuess<3>::
-  KA[GeometryInfo<3>::vertices_per_cell][3] =
-  {
-    {-0.250000, -0.250000, -0.250000},
-    { 0.250000, -0.250000, -0.250000},
-    {-0.250000,  0.250000, -0.250000},
-    { 0.250000,  0.250000, -0.250000},
-    {-0.250000, -0.250000,  0.250000},
-    { 0.250000, -0.250000,  0.250000},
-    {-0.250000,  0.250000,  0.250000},
-    { 0.250000,  0.250000,  0.250000}
-
-  };
-
-
-  template <>
-  const double
-  TransformR2UInitialGuess<3>::
-  Kb[GeometryInfo<3>::vertices_per_cell] =
-  {0.500000,0.250000,0.250000,0.000000,0.250000,0.000000,0.000000,-0.250000};
-
-  template<int dim, int spacedim>
-  Point<dim>
-  transform_real_to_unit_cell_initial_guess (const std::vector<Point<spacedim> > &vertex,
-                                             const Point<spacedim>               &p)
-  {
-    Point<dim> p_unit;
-
-    FullMatrix<double>  KA(GeometryInfo<dim>::vertices_per_cell, dim);
-    Vector <double>  Kb(GeometryInfo<dim>::vertices_per_cell);
-
-    KA.fill( (double *)(TransformR2UInitialGuess<dim>::KA) );
-    for (unsigned int i=0; i<GeometryInfo<dim>::vertices_per_cell; ++i)
-      Kb(i)=(TransformR2UInitialGuess<dim>::Kb)[i];
-
-    FullMatrix<double> Y(spacedim, GeometryInfo<dim>::vertices_per_cell);
-    for (unsigned int v=0; v<GeometryInfo<dim>::vertices_per_cell; v++)
-      for (unsigned int i=0; i<spacedim; ++i)
-        Y(i,v) = vertex[v][i];
-
-    FullMatrix<double> A(spacedim,dim);
-    Y.mmult(A,KA); // A = Y*KA
-    Vector< double > b(spacedim);
-    Y.vmult(b,Kb); // b = Y*Kb
-
-    for (unsigned int i=0; i<spacedim; ++i)
-      b(i) -= p[i];
-    b*=-1;
-
-    Vector< double > dest(dim);
-
-    FullMatrix<double> A_1(dim,spacedim);
-    if (dim<spacedim)
-      A_1.left_invert(A);
-    else
-      A_1.invert(A);
-
-    A_1.vmult(dest,b); //A^{-1}*b
-
-    for (unsigned int i=0; i<dim; ++i)
-      p_unit[i]=dest(i);
-
-    return p_unit;
-  }
-}
-
-
-
-
-template<int dim, int spacedim>
-Point<dim>
-MappingQ1<dim,spacedim>::
-transform_real_to_unit_cell (const typename Triangulation<dim,spacedim>::cell_iterator &cell,
-                             const Point<spacedim>                            &p) const
-{
-  // Use the exact formula if available
-  if (dim == spacedim && (dim == 1 || dim == 2))
-    {
-      // The dimension-dependent algorithms are much faster (about 25-45x in
-      // 2D) but fail most of the time when the given point (p) is not in the
-      // cell. The dimension-independent Newton algorithm given below is
-      // slower, but more robust (though it still sometimes fails). Therefore
-      // this function implements the following strategy based on the
-      // p's dimension:
-      //
-      // * In 1D this mapping is linear, so the mapping is always invertible
-      //   (and the exact formula is known) as long as the cell has non-zero
-      //   length.
-      // * In 2D the exact (quadratic) formula is called first. If either the
-      //   exact formula does not succeed (negative discriminant in the
-      //   quadratic formula) or succeeds but finds a solution outside of the
-      //   unit cell, then the Newton solver is called. The rationale for the
-      //   second choice is that the exact formula may provide two different
-      //   answers when mapping a point outside of the real cell, but the
-      //   Newton solver (if it converges) will only return one answer.
-      //   Otherwise the exact formula successfully found a point in the unit
-      //   cell and that value is returned.
-      // * In 3D there is no (known to the authors) exact formula, so the Newton
-      //   algorithm is used.
-      const std_cxx11::array<Point<spacedim>, GeometryInfo<dim>::vertices_per_cell>
-      vertices = this->get_vertices(cell);
-      try
-        {
-          Point<dim> point = internal::MappingQ1::transform_real_to_unit_cell(vertices, p);
-
-          if (dim == 1)
-            {
-              // formula not subject to any issues
-              return point;
-            }
-          else if (dim == 2)
-            {
-              // formula not guaranteed to work for points outside of the cell
-              const double eps = 1e-15;
-              if (-eps <= point(1) && point(1) <= 1 + eps &&
-                  -eps <= point(0) && point(0) <= 1 + eps)
-                {
-                  return point;
-                }
-            }
-          else
-            {
-              Assert(false, ExcInternalError());
-            }
-        }
-      catch (const typename Mapping<spacedim,spacedim>::ExcTransformationFailed &)
-        {
-          // continue on to the standard Newton code
-        }
-    }
-
-  // Find the initial value for the Newton iteration by a normal
-  // projection to the least square plane determined by the vertices
-  // of the cell
-  std::vector<Point<spacedim> > a;
-  this->compute_mapping_support_points (cell,a);
-  Assert(a.size() == GeometryInfo<dim>::vertices_per_cell,
-         ExcInternalError());
-
-  // if dim==1 there is nothing else to do to the initial value, and
-  // it is the answer
-  if (dim == 1)
-    return transform_real_to_unit_cell_initial_guess<dim,spacedim>(a,p);
-  else
-    {
-      const Point<dim> initial_p_unit =
-        transform_real_to_unit_cell_initial_guess<dim,spacedim>(a,p);
-
-      // 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 either the function computing an initial
-      // guess gave us a poor initial guess or for the following
-      // reason: we call this function here in the Q1 mapping to
-      // produce an initial guess for a higher order mapping, but 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
-
-      // perform the Newton iteration and
-      // return the result. note that this
-      // statement may throw an exception, which
-      // we simply pass up to the caller
-      return this->transform_real_to_unit_cell_internal(cell, p, initial_p_unit);
-    }
-}
-
-
-
-
-
 template<int dim, int spacedim>
 MappingQ1<dim,spacedim> *
 MappingQ1<dim,spacedim>::clone () const
index cefa8671c8fc191d5252430d87f8ec062ee935b0..034f2b4edb7741126d88bcc12dada25308e4b5aa 100644 (file)
@@ -1477,6 +1477,404 @@ transform_real_to_unit_cell_internal
 }
 
 
+namespace internal
+{
+  namespace MappingQ1
+  {
+    namespace
+    {
+
+      // These are left as templates on the spatial dimension (even though dim
+      // == spacedim must be true for them to make sense) because templates are
+      // expanded before the compiler eliminates code due to the 'if (dim ==
+      // spacedim)' statement (see the body of the general
+      // transform_real_to_unit_cell).
+      template<int spacedim>
+      Point<1>
+      transform_real_to_unit_cell
+      (const std_cxx11::array<Point<spacedim>, GeometryInfo<1>::vertices_per_cell> &vertices,
+       const Point<spacedim> &p)
+      {
+        Assert(spacedim == 1, ExcInternalError());
+        return Point<1>((p[0] - vertices[0](0))/(vertices[1](0) - vertices[0](0)));
+      }
+
+
+
+      template<int spacedim>
+      Point<2>
+      transform_real_to_unit_cell
+      (const std_cxx11::array<Point<spacedim>, GeometryInfo<2>::vertices_per_cell> &vertices,
+       const Point<spacedim> &p)
+      {
+        Assert(spacedim == 2, ExcInternalError());
+        const double x = p(0);
+        const double y = p(1);
+
+        const double x0 = vertices[0](0);
+        const double x1 = vertices[1](0);
+        const double x2 = vertices[2](0);
+        const double x3 = vertices[3](0);
+
+        const double y0 = vertices[0](1);
+        const double y1 = vertices[1](1);
+        const double y2 = vertices[2](1);
+        const double y3 = vertices[3](1);
+
+        const double a = (x1 - x3)*(y0 - y2) - (x0 - x2)*(y1 - y3);
+        const double b = -(x0 - x1 - x2 + x3)*y + (x - 2*x1 + x3)*y0 - (x - 2*x0 + x2)*y1
+                         - (x - x1)*y2 + (x - x0)*y3;
+        const double c = (x0 - x1)*y - (x - x1)*y0 + (x - x0)*y1;
+
+        const double discriminant = b*b - 4*a*c;
+        // exit if the point is not in the cell (this is the only case where the
+        // discriminant is negative)
+        if (discriminant < 0.0)
+          {
+            AssertThrow (false,
+                         (typename Mapping<spacedim,spacedim>::ExcTransformationFailed()));
+          }
+
+        double eta1;
+        double eta2;
+        // special case #1: if a is zero, then use the linear formula
+        if (a == 0.0 && b != 0.0)
+          {
+            eta1 = -c/b;
+            eta2 = -c/b;
+          }
+        // special case #2: if c is very small:
+        else if (std::abs(c/b) < 1e-12)
+          {
+            eta1 = (-b - std::sqrt(discriminant)) / (2*a);
+            eta2 = (-b + std::sqrt(discriminant)) / (2*a);
+          }
+        // finally, use the numerically stable version of the quadratic formula:
+        else
+          {
+            eta1 = 2*c / (-b - std::sqrt(discriminant));
+            eta2 = 2*c / (-b + std::sqrt(discriminant));
+          }
+        // pick the one closer to the center of the cell.
+        const double eta = (std::abs(eta1 - 0.5) < std::abs(eta2 - 0.5)) ? eta1 : eta2;
+
+        /*
+         * There are two ways to compute xi from eta, but either one may have a
+         * zero denominator.
+         */
+        const double subexpr0 = -eta*x2 + x0*(eta - 1);
+        const double xi_denominator0 = eta*x3 - x1*(eta - 1) + subexpr0;
+        const double max_x = std::max(std::max(std::abs(x0), std::abs(x1)),
+                                      std::max(std::abs(x2), std::abs(x3)));
+
+        if (std::abs(xi_denominator0) > 1e-10*max_x)
+          {
+            const double xi = (x + subexpr0)/xi_denominator0;
+            return Point<2>(xi, eta);
+          }
+        else
+          {
+            const double max_y = std::max(std::max(std::abs(y0), std::abs(y1)),
+                                          std::max(std::abs(y2), std::abs(y3)));
+            const double subexpr1 = -eta*y2 + y0*(eta - 1);
+            const double xi_denominator1 = eta*y3 - y1*(eta - 1) + subexpr1;
+            if (std::abs(xi_denominator1) > 1e-10*max_y)
+              {
+                const double xi = (subexpr1 + y)/xi_denominator1;
+                return Point<2>(xi, eta);
+              }
+            else // give up and try Newton iteration
+              {
+                AssertThrow (false,
+                             (typename Mapping<spacedim,spacedim>::ExcTransformationFailed()));
+              }
+          }
+        // bogus return to placate compiler. It should not be possible to get
+        // here.
+        Assert(false, ExcInternalError());
+        return Point<2>(std::numeric_limits<double>::quiet_NaN(),
+                        std::numeric_limits<double>::quiet_NaN());
+      }
+
+
+
+      template<int spacedim>
+      Point<3>
+      transform_real_to_unit_cell
+      (const std_cxx11::array<Point<spacedim>, GeometryInfo<3>::vertices_per_cell> &/*vertices*/,
+       const Point<spacedim> &/*p*/)
+      {
+        // It should not be possible to get here
+        Assert(false, ExcInternalError());
+        return Point<3>();
+      }
+
+
+
+      /**
+       * Compute an initial guess to pass to the Newton method in
+       * transform_real_to_unit_cell.  For the initial guess we proceed in the
+       * following way:
+       * <ul>
+       * <li> find the least square dim-dimensional plane approximating the cell
+       * vertices, i.e. we find an affine map A x_hat + b from the reference cell
+       * to the real space.
+       * <li> Solve the equation A x_hat + b = p for x_hat
+       * <li> This x_hat is the initial solution used for the Newton Method.
+       * </ul>
+       *
+       * @note if dim<spacedim we first project p onto the plane.
+       *
+       * @note if dim==1 (for any spacedim) the initial guess is the exact
+       * solution and no Newton iteration is needed.
+       *
+       * Some details about how we compute the least square plane. We look
+       * for a spacedim x (dim + 1) matrix X such that X * M = Y where M is
+       * a (dim+1) x n_vertices matrix and Y a spacedim x n_vertices.  And:
+       * The i-th column of M is unit_vertex[i] and the last row all
+       * 1's. The i-th column of Y is real_vertex[i].  If we split X=[A|b],
+       * the least square approx is A x_hat+b Classically X = Y * (M^t (M
+       * M^t)^{-1}) Let K = M^t * (M M^t)^{-1} = [KA Kb] this can be
+       * precomputed, and that is exactly what we do.  Finally A = Y*KA and
+       * b = Y*Kb.
+       */
+      template <int dim>
+      struct TransformR2UInitialGuess
+      {
+        static const double KA[GeometryInfo<dim>::vertices_per_cell][dim];
+        static const double Kb[GeometryInfo<dim>::vertices_per_cell];
+      };
+
+
+      /*
+        Octave code:
+        M=[0 1; 1 1];
+        K1 = transpose(M) * inverse (M*transpose(M));
+        printf ("{%f, %f},\n", K1' );
+      */
+      template <>
+      const double
+      TransformR2UInitialGuess<1>::
+      KA[GeometryInfo<1>::vertices_per_cell][1] =
+      {
+        {-1.000000},
+        {1.000000}
+      };
+
+      template <>
+      const double
+      TransformR2UInitialGuess<1>::
+      Kb[GeometryInfo<1>::vertices_per_cell] = {1.000000, 0.000000};
+
+
+      /*
+        Octave code:
+        M=[0 1 0 1;0 0 1 1;1 1 1 1];
+        K2 = transpose(M) * inverse (M*transpose(M));
+        printf ("{%f, %f, %f},\n", K2' );
+      */
+      template <>
+      const double
+      TransformR2UInitialGuess<2>::
+      KA[GeometryInfo<2>::vertices_per_cell][2] =
+      {
+        {-0.500000, -0.500000},
+        { 0.500000, -0.500000},
+        {-0.500000,  0.500000},
+        { 0.500000,  0.500000}
+      };
+
+      /*
+        Octave code:
+        M=[0 1 0 1 0 1 0 1;0 0 1 1 0 0 1 1; 0 0 0 0 1 1 1 1; 1 1 1 1 1 1 1 1];
+        K3 = transpose(M) * inverse (M*transpose(M))
+        printf ("{%f, %f, %f, %f},\n", K3' );
+      */
+      template <>
+      const double
+      TransformR2UInitialGuess<2>::
+      Kb[GeometryInfo<2>::vertices_per_cell] =
+      {0.750000,0.250000,0.250000,-0.250000 };
+
+
+      template <>
+      const double
+      TransformR2UInitialGuess<3>::
+      KA[GeometryInfo<3>::vertices_per_cell][3] =
+      {
+        {-0.250000, -0.250000, -0.250000},
+        { 0.250000, -0.250000, -0.250000},
+        {-0.250000,  0.250000, -0.250000},
+        { 0.250000,  0.250000, -0.250000},
+        {-0.250000, -0.250000,  0.250000},
+        { 0.250000, -0.250000,  0.250000},
+        {-0.250000,  0.250000,  0.250000},
+        { 0.250000,  0.250000,  0.250000}
+
+      };
+
+
+      template <>
+      const double
+      TransformR2UInitialGuess<3>::
+      Kb[GeometryInfo<3>::vertices_per_cell] =
+      {0.500000,0.250000,0.250000,0.000000,0.250000,0.000000,0.000000,-0.250000};
+
+      template<int dim, int spacedim>
+      Point<dim>
+      transform_real_to_unit_cell_initial_guess (const std::vector<Point<spacedim> > &vertex,
+                                                 const Point<spacedim>               &p)
+      {
+        Point<dim> p_unit;
+
+        dealii::FullMatrix<double>  KA(GeometryInfo<dim>::vertices_per_cell, dim);
+        dealii::Vector <double>  Kb(GeometryInfo<dim>::vertices_per_cell);
+
+        KA.fill( (double *)(TransformR2UInitialGuess<dim>::KA) );
+        for (unsigned int i=0; i<GeometryInfo<dim>::vertices_per_cell; ++i)
+          Kb(i) = TransformR2UInitialGuess<dim>::Kb[i];
+
+        FullMatrix<double> Y(spacedim, GeometryInfo<dim>::vertices_per_cell);
+        for (unsigned int v=0; v<GeometryInfo<dim>::vertices_per_cell; v++)
+          for (unsigned int i=0; i<spacedim; ++i)
+            Y(i,v) = vertex[v][i];
+
+        FullMatrix<double> A(spacedim,dim);
+        Y.mmult(A,KA); // A = Y*KA
+        dealii::Vector<double> b(spacedim);
+        Y.vmult(b,Kb); // b = Y*Kb
+
+        for (unsigned int i=0; i<spacedim; ++i)
+          b(i) -= p[i];
+        b*=-1;
+
+        dealii::Vector<double> dest(dim);
+
+        FullMatrix<double> A_1(dim,spacedim);
+        if (dim<spacedim)
+          A_1.left_invert(A);
+        else
+          A_1.invert(A);
+
+        A_1.vmult(dest,b); //A^{-1}*b
+
+        for (unsigned int i=0; i<dim; ++i)
+          p_unit[i]=dest(i);
+
+        return p_unit;
+      }
+    }
+  }
+}
+
+
+
+
+template<int dim, int spacedim>
+Point<dim>
+MappingQGeneric<dim,spacedim>::
+transform_real_to_unit_cell (const typename Triangulation<dim,spacedim>::cell_iterator &cell,
+                             const Point<spacedim>                            &p) const
+{
+  // Use an exact formula if one is available. this is only the case
+  // for Q1 mappings in 1d, and in 2d if dim==spacedim
+  if ((polynomial_degree == 1) &&
+      ((dim == 1)
+       ||
+       ((dim == 2) && (dim == spacedim))))
+    {
+      // The dimension-dependent algorithms are much faster (about 25-45x in
+      // 2D) but fail most of the time when the given point (p) is not in the
+      // cell. The dimension-independent Newton algorithm given below is
+      // slower, but more robust (though it still sometimes fails). Therefore
+      // this function implements the following strategy based on the
+      // p's dimension:
+      //
+      // * In 1D this mapping is linear, so the mapping is always invertible
+      //   (and the exact formula is known) as long as the cell has non-zero
+      //   length.
+      // * In 2D the exact (quadratic) formula is called first. If either the
+      //   exact formula does not succeed (negative discriminant in the
+      //   quadratic formula) or succeeds but finds a solution outside of the
+      //   unit cell, then the Newton solver is called. The rationale for the
+      //   second choice is that the exact formula may provide two different
+      //   answers when mapping a point outside of the real cell, but the
+      //   Newton solver (if it converges) will only return one answer.
+      //   Otherwise the exact formula successfully found a point in the unit
+      //   cell and that value is returned.
+      // * In 3D there is no (known to the authors) exact formula, so the Newton
+      //   algorithm is used.
+      const std_cxx11::array<Point<spacedim>, GeometryInfo<dim>::vertices_per_cell>
+      vertices = this->get_vertices(cell);
+      try
+        {
+          switch (dim)
+            {
+            case 1:
+            {
+              // formula not subject to any issues in 1d
+              if (spacedim == 1)
+                return internal::MappingQ1::transform_real_to_unit_cell(vertices, p);
+              else
+                {
+                  const std::vector<Point<spacedim> > a (vertices.begin(),
+                                                         vertices.end());
+                  return internal::MappingQ1::transform_real_to_unit_cell_initial_guess<dim,spacedim>(a,p);
+                }
+            }
+
+            case 2:
+            {
+              const Point<dim> point
+                = internal::MappingQ1::transform_real_to_unit_cell(vertices, p);
+
+              // formula not guaranteed to work for points outside of
+              // the cell. only take the computed point if it lies
+              // inside the reference cell
+              const double eps = 1e-15;
+              if (-eps <= point(1) && point(1) <= 1 + eps &&
+                  -eps <= point(0) && point(0) <= 1 + eps)
+                {
+                  return point;
+                }
+              else
+                break;
+            }
+
+            default:
+            {
+              // we should get here, based on the if-condition at the top
+              Assert(false, ExcInternalError());
+            }
+            }
+        }
+      catch (const typename Mapping<spacedim,spacedim>::ExcTransformationFailed &)
+        {
+          // simply fall through and continue on to the standard Newton code
+        }
+    }
+  else
+    {
+      // we can't use an explicit formula,
+    }
+
+
+  // Find the initial value for the Newton iteration by a normal
+  // projection to the least square plane determined by the vertices
+  // of the cell
+  std::vector<Point<spacedim> > a;
+  compute_mapping_support_points (cell,a);
+  Assert(a.size() == GeometryInfo<dim>::vertices_per_cell,
+         ExcInternalError());
+  const Point<dim> initial_p_unit =
+    internal::MappingQ1::transform_real_to_unit_cell_initial_guess<dim,spacedim>(a,p);
+
+  // perform the Newton iteration and return the result. note that
+  // this statement may throw an exception, which we simply pass up to
+  // the caller
+  return this->transform_real_to_unit_cell_internal(cell, p, initial_p_unit);
+}
+
 
 
 template<int dim, int spacedim>

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.