]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Remove the internal::MappingQGeneric namespace. 1764/head
authorDavid Wells <wellsd2@rpi.edu>
Tue, 20 Oct 2015 15:15:08 +0000 (11:15 -0400)
committerDavid Wells <wellsd2@rpi.edu>
Tue, 20 Oct 2015 15:33:32 +0000 (11:33 -0400)
The diff git came up with is much more complex than what was actually
changed. This commit moves all the functions that were in
internal::MappingQGeneric into internal::MappingQ1 and places the new
combined namespace at the top of the file.

source/fe/mapping_q_generic.cc

index 91db94cfc0bbd2fb9fe7428659a5901fdfd0f35c..7095eba7e982b3904def498ffcd57d1977cd302d 100644 (file)
 
 DEAL_II_NAMESPACE_OPEN
 
+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 spacedim>
+      void
+      compute_shape_function_values (const unsigned int            n_shape_functions,
+                                     const std::vector<Point<1> > &unit_points,
+                                     typename dealii::MappingQGeneric<1,spacedim>::InternalData &data)
+      {
+        (void)n_shape_functions;
+        const unsigned int n_points=unit_points.size();
+        for (unsigned int k = 0 ; k < n_points ; ++k)
+          {
+            double x = unit_points[k](0);
+
+            if (data.shape_values.size()!=0)
+              {
+                Assert(data.shape_values.size()==n_shape_functions*n_points,
+                       ExcInternalError());
+                data.shape(k,0) = 1.-x;
+                data.shape(k,1) = x;
+              }
+            if (data.shape_derivatives.size()!=0)
+              {
+                Assert(data.shape_derivatives.size()==n_shape_functions*n_points,
+                       ExcInternalError());
+                data.derivative(k,0)[0] = -1.;
+                data.derivative(k,1)[0] = 1.;
+              }
+            if (data.shape_second_derivatives.size()!=0)
+              {
+                // the following may or may not
+                // work if dim != spacedim
+                Assert (spacedim == 1, ExcNotImplemented());
+
+                Assert(data.shape_second_derivatives.size()==n_shape_functions*n_points,
+                       ExcInternalError());
+                data.second_derivative(k,0)[0][0] = 0;
+                data.second_derivative(k,1)[0][0] = 0;
+              }
+            if (data.shape_third_derivatives.size()!=0)
+              {
+                // if lower order derivative don't work, neither should this
+                Assert (spacedim == 1, ExcNotImplemented());
+
+                Assert(data.shape_third_derivatives.size()==n_shape_functions*n_points,
+                       ExcInternalError());
+
+                Tensor<3,1> zero;
+                data.third_derivative(k,0) = zero;
+                data.third_derivative(k,1) = zero;
+              }
+            if (data.shape_fourth_derivatives.size()!=0)
+              {
+                // if lower order derivative don't work, neither should this
+                Assert (spacedim == 1, ExcNotImplemented());
+
+                Assert(data.shape_fourth_derivatives.size()==n_shape_functions*n_points,
+                       ExcInternalError());
+
+                Tensor<4,1> zero;
+                data.fourth_derivative(k,0) = zero;
+                data.fourth_derivative(k,1) = zero;
+              }
+          }
+      }
+
+
+      template <int spacedim>
+      void
+      compute_shape_function_values (const unsigned int            n_shape_functions,
+                                     const std::vector<Point<2> > &unit_points,
+                                     typename dealii::MappingQGeneric<2,spacedim>::InternalData &data)
+      {
+        (void)n_shape_functions;
+        const unsigned int n_points=unit_points.size();
+        for (unsigned int k = 0 ; k < n_points ; ++k)
+          {
+            double x = unit_points[k](0);
+            double y = unit_points[k](1);
+
+            if (data.shape_values.size()!=0)
+              {
+                Assert(data.shape_values.size()==n_shape_functions*n_points,
+                       ExcInternalError());
+                data.shape(k,0) = (1.-x)*(1.-y);
+                data.shape(k,1) = x*(1.-y);
+                data.shape(k,2) = (1.-x)*y;
+                data.shape(k,3) = x*y;
+              }
+            if (data.shape_derivatives.size()!=0)
+              {
+                Assert(data.shape_derivatives.size()==n_shape_functions*n_points,
+                       ExcInternalError());
+                data.derivative(k,0)[0] = (y-1.);
+                data.derivative(k,1)[0] = (1.-y);
+                data.derivative(k,2)[0] = -y;
+                data.derivative(k,3)[0] = y;
+                data.derivative(k,0)[1] = (x-1.);
+                data.derivative(k,1)[1] = -x;
+                data.derivative(k,2)[1] = (1.-x);
+                data.derivative(k,3)[1] = x;
+              }
+            if (data.shape_second_derivatives.size()!=0)
+              {
+                Assert(data.shape_second_derivatives.size()==n_shape_functions*n_points,
+                       ExcInternalError());
+                data.second_derivative(k,0)[0][0] = 0;
+                data.second_derivative(k,1)[0][0] = 0;
+                data.second_derivative(k,2)[0][0] = 0;
+                data.second_derivative(k,3)[0][0] = 0;
+                data.second_derivative(k,0)[0][1] = 1.;
+                data.second_derivative(k,1)[0][1] = -1.;
+                data.second_derivative(k,2)[0][1] = -1.;
+                data.second_derivative(k,3)[0][1] = 1.;
+                data.second_derivative(k,0)[1][0] = 1.;
+                data.second_derivative(k,1)[1][0] = -1.;
+                data.second_derivative(k,2)[1][0] = -1.;
+                data.second_derivative(k,3)[1][0] = 1.;
+                data.second_derivative(k,0)[1][1] = 0;
+                data.second_derivative(k,1)[1][1] = 0;
+                data.second_derivative(k,2)[1][1] = 0;
+                data.second_derivative(k,3)[1][1] = 0;
+              }
+            if (data.shape_third_derivatives.size()!=0)
+              {
+                Assert(data.shape_third_derivatives.size()==n_shape_functions*n_points,
+                       ExcInternalError());
+
+                Tensor<3,2> zero;
+                for (unsigned int i=0; i<4; ++i)
+                  data.third_derivative(k,i) = zero;
+              }
+            if (data.shape_fourth_derivatives.size()!=0)
+              {
+                Assert(data.shape_fourth_derivatives.size()==n_shape_functions*n_points,
+                       ExcInternalError());
+                Tensor<4,2> zero;
+                for (unsigned int i=0; i<4; ++i)
+                  data.fourth_derivative(k,i) = zero;
+              }
+          }
+      }
+
+
+
+      template <int spacedim>
+      void
+      compute_shape_function_values (const unsigned int            n_shape_functions,
+                                     const std::vector<Point<3> > &unit_points,
+                                     typename dealii::MappingQGeneric<3,spacedim>::InternalData &data)
+      {
+        (void)n_shape_functions;
+        const unsigned int n_points=unit_points.size();
+        for (unsigned int k = 0 ; k < n_points ; ++k)
+          {
+            double x = unit_points[k](0);
+            double y = unit_points[k](1);
+            double z = unit_points[k](2);
+
+            if (data.shape_values.size()!=0)
+              {
+                Assert(data.shape_values.size()==n_shape_functions*n_points,
+                       ExcInternalError());
+                data.shape(k,0) = (1.-x)*(1.-y)*(1.-z);
+                data.shape(k,1) = x*(1.-y)*(1.-z);
+                data.shape(k,2) = (1.-x)*y*(1.-z);
+                data.shape(k,3) = x*y*(1.-z);
+                data.shape(k,4) = (1.-x)*(1.-y)*z;
+                data.shape(k,5) = x*(1.-y)*z;
+                data.shape(k,6) = (1.-x)*y*z;
+                data.shape(k,7) = x*y*z;
+              }
+            if (data.shape_derivatives.size()!=0)
+              {
+                Assert(data.shape_derivatives.size()==n_shape_functions*n_points,
+                       ExcInternalError());
+                data.derivative(k,0)[0] = (y-1.)*(1.-z);
+                data.derivative(k,1)[0] = (1.-y)*(1.-z);
+                data.derivative(k,2)[0] = -y*(1.-z);
+                data.derivative(k,3)[0] = y*(1.-z);
+                data.derivative(k,4)[0] = (y-1.)*z;
+                data.derivative(k,5)[0] = (1.-y)*z;
+                data.derivative(k,6)[0] = -y*z;
+                data.derivative(k,7)[0] = y*z;
+                data.derivative(k,0)[1] = (x-1.)*(1.-z);
+                data.derivative(k,1)[1] = -x*(1.-z);
+                data.derivative(k,2)[1] = (1.-x)*(1.-z);
+                data.derivative(k,3)[1] = x*(1.-z);
+                data.derivative(k,4)[1] = (x-1.)*z;
+                data.derivative(k,5)[1] = -x*z;
+                data.derivative(k,6)[1] = (1.-x)*z;
+                data.derivative(k,7)[1] = x*z;
+                data.derivative(k,0)[2] = (x-1)*(1.-y);
+                data.derivative(k,1)[2] = x*(y-1.);
+                data.derivative(k,2)[2] = (x-1.)*y;
+                data.derivative(k,3)[2] = -x*y;
+                data.derivative(k,4)[2] = (1.-x)*(1.-y);
+                data.derivative(k,5)[2] = x*(1.-y);
+                data.derivative(k,6)[2] = (1.-x)*y;
+                data.derivative(k,7)[2] = x*y;
+              }
+            if (data.shape_second_derivatives.size()!=0)
+              {
+                // the following may or may not
+                // work if dim != spacedim
+                Assert (spacedim == 3, ExcNotImplemented());
+
+                Assert(data.shape_second_derivatives.size()==n_shape_functions*n_points,
+                       ExcInternalError());
+                data.second_derivative(k,0)[0][0] = 0;
+                data.second_derivative(k,1)[0][0] = 0;
+                data.second_derivative(k,2)[0][0] = 0;
+                data.second_derivative(k,3)[0][0] = 0;
+                data.second_derivative(k,4)[0][0] = 0;
+                data.second_derivative(k,5)[0][0] = 0;
+                data.second_derivative(k,6)[0][0] = 0;
+                data.second_derivative(k,7)[0][0] = 0;
+                data.second_derivative(k,0)[1][1] = 0;
+                data.second_derivative(k,1)[1][1] = 0;
+                data.second_derivative(k,2)[1][1] = 0;
+                data.second_derivative(k,3)[1][1] = 0;
+                data.second_derivative(k,4)[1][1] = 0;
+                data.second_derivative(k,5)[1][1] = 0;
+                data.second_derivative(k,6)[1][1] = 0;
+                data.second_derivative(k,7)[1][1] = 0;
+                data.second_derivative(k,0)[2][2] = 0;
+                data.second_derivative(k,1)[2][2] = 0;
+                data.second_derivative(k,2)[2][2] = 0;
+                data.second_derivative(k,3)[2][2] = 0;
+                data.second_derivative(k,4)[2][2] = 0;
+                data.second_derivative(k,5)[2][2] = 0;
+                data.second_derivative(k,6)[2][2] = 0;
+                data.second_derivative(k,7)[2][2] = 0;
+
+                data.second_derivative(k,0)[0][1] = (1.-z);
+                data.second_derivative(k,1)[0][1] = -(1.-z);
+                data.second_derivative(k,2)[0][1] = -(1.-z);
+                data.second_derivative(k,3)[0][1] = (1.-z);
+                data.second_derivative(k,4)[0][1] = z;
+                data.second_derivative(k,5)[0][1] = -z;
+                data.second_derivative(k,6)[0][1] = -z;
+                data.second_derivative(k,7)[0][1] = z;
+                data.second_derivative(k,0)[1][0] = (1.-z);
+                data.second_derivative(k,1)[1][0] = -(1.-z);
+                data.second_derivative(k,2)[1][0] = -(1.-z);
+                data.second_derivative(k,3)[1][0] = (1.-z);
+                data.second_derivative(k,4)[1][0] = z;
+                data.second_derivative(k,5)[1][0] = -z;
+                data.second_derivative(k,6)[1][0] = -z;
+                data.second_derivative(k,7)[1][0] = z;
+
+                data.second_derivative(k,0)[0][2] = (1.-y);
+                data.second_derivative(k,1)[0][2] = -(1.-y);
+                data.second_derivative(k,2)[0][2] = y;
+                data.second_derivative(k,3)[0][2] = -y;
+                data.second_derivative(k,4)[0][2] = -(1.-y);
+                data.second_derivative(k,5)[0][2] = (1.-y);
+                data.second_derivative(k,6)[0][2] = -y;
+                data.second_derivative(k,7)[0][2] = y;
+                data.second_derivative(k,0)[2][0] = (1.-y);
+                data.second_derivative(k,1)[2][0] = -(1.-y);
+                data.second_derivative(k,2)[2][0] = y;
+                data.second_derivative(k,3)[2][0] = -y;
+                data.second_derivative(k,4)[2][0] = -(1.-y);
+                data.second_derivative(k,5)[2][0] = (1.-y);
+                data.second_derivative(k,6)[2][0] = -y;
+                data.second_derivative(k,7)[2][0] = y;
+
+                data.second_derivative(k,0)[1][2] = (1.-x);
+                data.second_derivative(k,1)[1][2] = x;
+                data.second_derivative(k,2)[1][2] = -(1.-x);
+                data.second_derivative(k,3)[1][2] = -x;
+                data.second_derivative(k,4)[1][2] = -(1.-x);
+                data.second_derivative(k,5)[1][2] = -x;
+                data.second_derivative(k,6)[1][2] = (1.-x);
+                data.second_derivative(k,7)[1][2] = x;
+                data.second_derivative(k,0)[2][1] = (1.-x);
+                data.second_derivative(k,1)[2][1] = x;
+                data.second_derivative(k,2)[2][1] = -(1.-x);
+                data.second_derivative(k,3)[2][1] = -x;
+                data.second_derivative(k,4)[2][1] = -(1.-x);
+                data.second_derivative(k,5)[2][1] = -x;
+                data.second_derivative(k,6)[2][1] = (1.-x);
+                data.second_derivative(k,7)[2][1] = x;
+              }
+            if (data.shape_third_derivatives.size()!=0)
+              {
+                // if lower order derivative don't work, neither should this
+                Assert (spacedim == 3, ExcNotImplemented());
+
+                Assert(data.shape_third_derivatives.size()==n_shape_functions*n_points,
+                       ExcInternalError());
+
+                for (unsigned int i=0; i<3; ++i)
+                  for (unsigned int j=0; j<3; ++j)
+                    for (unsigned int l=0; l<3; ++l)
+                      if ((i==j)||(j==l)||(l==i))
+                        {
+                          for (unsigned int m=0; m<8; ++m)
+                            data.third_derivative(k,m)[i][j][l] = 0;
+                        }
+                      else
+                        {
+                          data.third_derivative(k,0)[i][j][l] = -1.;
+                          data.third_derivative(k,1)[i][j][l] = 1.;
+                          data.third_derivative(k,2)[i][j][l] = 1.;
+                          data.third_derivative(k,3)[i][j][l] = -1.;
+                          data.third_derivative(k,4)[i][j][l] = 1.;
+                          data.third_derivative(k,5)[i][j][l] = -1.;
+                          data.third_derivative(k,6)[i][j][l] = -1.;
+                          data.third_derivative(k,7)[i][j][l] = 1.;
+                        }
+
+              }
+            if (data.shape_fourth_derivatives.size()!=0)
+              {
+                // if lower order derivative don't work, neither should this
+                Assert (spacedim == 3, ExcNotImplemented());
+
+                Assert(data.shape_fourth_derivatives.size()==n_shape_functions*n_points,
+                       ExcInternalError());
+                Tensor<4,3> zero;
+                for (unsigned int i=0; i<8; ++i)
+                  data.fourth_derivative(k,i) = zero;
+              }
+          }
+      }
+    }
+  }
+}
+
+
+
+
 
 template<int dim, int spacedim>
 MappingQGeneric<dim,spacedim>::InternalData::InternalData (const unsigned int polynomial_degree)
@@ -188,348 +806,15 @@ initialize_face (const UpdateFlags      update_flags,
                   std::fill (unit_tangentials[i].begin(),
                              unit_tangentials[i].end(), tang1);
                   std::fill (unit_tangentials[nfaces+i].begin(),
-                             unit_tangentials[nfaces+i].end(), tang2);
-                }
-            }
-        }
-    }
-}
-
-
-
-namespace internal
-{
-  namespace MappingQGeneric
-  {
-    template <int spacedim>
-    void
-    compute_shape_function_values (const unsigned int            n_shape_functions,
-                                   const std::vector<Point<1> > &unit_points,
-                                   typename dealii::MappingQGeneric<1,spacedim>::InternalData &data)
-    {
-      (void)n_shape_functions;
-      const unsigned int n_points=unit_points.size();
-      for (unsigned int k = 0 ; k < n_points ; ++k)
-        {
-          double x = unit_points[k](0);
-
-          if (data.shape_values.size()!=0)
-            {
-              Assert(data.shape_values.size()==n_shape_functions*n_points,
-                     ExcInternalError());
-              data.shape(k,0) = 1.-x;
-              data.shape(k,1) = x;
-            }
-          if (data.shape_derivatives.size()!=0)
-            {
-              Assert(data.shape_derivatives.size()==n_shape_functions*n_points,
-                     ExcInternalError());
-              data.derivative(k,0)[0] = -1.;
-              data.derivative(k,1)[0] = 1.;
-            }
-          if (data.shape_second_derivatives.size()!=0)
-            {
-              // the following may or may not
-              // work if dim != spacedim
-              Assert (spacedim == 1, ExcNotImplemented());
-
-              Assert(data.shape_second_derivatives.size()==n_shape_functions*n_points,
-                     ExcInternalError());
-              data.second_derivative(k,0)[0][0] = 0;
-              data.second_derivative(k,1)[0][0] = 0;
-            }
-          if (data.shape_third_derivatives.size()!=0)
-            {
-              // if lower order derivative don't work, neither should this
-              Assert (spacedim == 1, ExcNotImplemented());
-
-              Assert(data.shape_third_derivatives.size()==n_shape_functions*n_points,
-                     ExcInternalError());
-
-              Tensor<3,1> zero;
-              data.third_derivative(k,0) = zero;
-              data.third_derivative(k,1) = zero;
-            }
-          if (data.shape_fourth_derivatives.size()!=0)
-            {
-              // if lower order derivative don't work, neither should this
-              Assert (spacedim == 1, ExcNotImplemented());
-
-              Assert(data.shape_fourth_derivatives.size()==n_shape_functions*n_points,
-                     ExcInternalError());
-
-              Tensor<4,1> zero;
-              data.fourth_derivative(k,0) = zero;
-              data.fourth_derivative(k,1) = zero;
-            }
-        }
-    }
-
-
-    template <int spacedim>
-    void
-    compute_shape_function_values (const unsigned int            n_shape_functions,
-                                   const std::vector<Point<2> > &unit_points,
-                                   typename dealii::MappingQGeneric<2,spacedim>::InternalData &data)
-    {
-      (void)n_shape_functions;
-      const unsigned int n_points=unit_points.size();
-      for (unsigned int k = 0 ; k < n_points ; ++k)
-        {
-          double x = unit_points[k](0);
-          double y = unit_points[k](1);
-
-          if (data.shape_values.size()!=0)
-            {
-              Assert(data.shape_values.size()==n_shape_functions*n_points,
-                     ExcInternalError());
-              data.shape(k,0) = (1.-x)*(1.-y);
-              data.shape(k,1) = x*(1.-y);
-              data.shape(k,2) = (1.-x)*y;
-              data.shape(k,3) = x*y;
-            }
-          if (data.shape_derivatives.size()!=0)
-            {
-              Assert(data.shape_derivatives.size()==n_shape_functions*n_points,
-                     ExcInternalError());
-              data.derivative(k,0)[0] = (y-1.);
-              data.derivative(k,1)[0] = (1.-y);
-              data.derivative(k,2)[0] = -y;
-              data.derivative(k,3)[0] = y;
-              data.derivative(k,0)[1] = (x-1.);
-              data.derivative(k,1)[1] = -x;
-              data.derivative(k,2)[1] = (1.-x);
-              data.derivative(k,3)[1] = x;
-            }
-          if (data.shape_second_derivatives.size()!=0)
-            {
-              Assert(data.shape_second_derivatives.size()==n_shape_functions*n_points,
-                     ExcInternalError());
-              data.second_derivative(k,0)[0][0] = 0;
-              data.second_derivative(k,1)[0][0] = 0;
-              data.second_derivative(k,2)[0][0] = 0;
-              data.second_derivative(k,3)[0][0] = 0;
-              data.second_derivative(k,0)[0][1] = 1.;
-              data.second_derivative(k,1)[0][1] = -1.;
-              data.second_derivative(k,2)[0][1] = -1.;
-              data.second_derivative(k,3)[0][1] = 1.;
-              data.second_derivative(k,0)[1][0] = 1.;
-              data.second_derivative(k,1)[1][0] = -1.;
-              data.second_derivative(k,2)[1][0] = -1.;
-              data.second_derivative(k,3)[1][0] = 1.;
-              data.second_derivative(k,0)[1][1] = 0;
-              data.second_derivative(k,1)[1][1] = 0;
-              data.second_derivative(k,2)[1][1] = 0;
-              data.second_derivative(k,3)[1][1] = 0;
-            }
-          if (data.shape_third_derivatives.size()!=0)
-            {
-              Assert(data.shape_third_derivatives.size()==n_shape_functions*n_points,
-                     ExcInternalError());
-
-              Tensor<3,2> zero;
-              for (unsigned int i=0; i<4; ++i)
-                data.third_derivative(k,i) = zero;
-            }
-          if (data.shape_fourth_derivatives.size()!=0)
-            {
-              Assert(data.shape_fourth_derivatives.size()==n_shape_functions*n_points,
-                     ExcInternalError());
-              Tensor<4,2> zero;
-              for (unsigned int i=0; i<4; ++i)
-                data.fourth_derivative(k,i) = zero;
-            }
-        }
-    }
-
-
-
-    template <int spacedim>
-    void
-    compute_shape_function_values (const unsigned int            n_shape_functions,
-                                   const std::vector<Point<3> > &unit_points,
-                                   typename dealii::MappingQGeneric<3,spacedim>::InternalData &data)
-    {
-      (void)n_shape_functions;
-      const unsigned int n_points=unit_points.size();
-      for (unsigned int k = 0 ; k < n_points ; ++k)
-        {
-          double x = unit_points[k](0);
-          double y = unit_points[k](1);
-          double z = unit_points[k](2);
-
-          if (data.shape_values.size()!=0)
-            {
-              Assert(data.shape_values.size()==n_shape_functions*n_points,
-                     ExcInternalError());
-              data.shape(k,0) = (1.-x)*(1.-y)*(1.-z);
-              data.shape(k,1) = x*(1.-y)*(1.-z);
-              data.shape(k,2) = (1.-x)*y*(1.-z);
-              data.shape(k,3) = x*y*(1.-z);
-              data.shape(k,4) = (1.-x)*(1.-y)*z;
-              data.shape(k,5) = x*(1.-y)*z;
-              data.shape(k,6) = (1.-x)*y*z;
-              data.shape(k,7) = x*y*z;
-            }
-          if (data.shape_derivatives.size()!=0)
-            {
-              Assert(data.shape_derivatives.size()==n_shape_functions*n_points,
-                     ExcInternalError());
-              data.derivative(k,0)[0] = (y-1.)*(1.-z);
-              data.derivative(k,1)[0] = (1.-y)*(1.-z);
-              data.derivative(k,2)[0] = -y*(1.-z);
-              data.derivative(k,3)[0] = y*(1.-z);
-              data.derivative(k,4)[0] = (y-1.)*z;
-              data.derivative(k,5)[0] = (1.-y)*z;
-              data.derivative(k,6)[0] = -y*z;
-              data.derivative(k,7)[0] = y*z;
-              data.derivative(k,0)[1] = (x-1.)*(1.-z);
-              data.derivative(k,1)[1] = -x*(1.-z);
-              data.derivative(k,2)[1] = (1.-x)*(1.-z);
-              data.derivative(k,3)[1] = x*(1.-z);
-              data.derivative(k,4)[1] = (x-1.)*z;
-              data.derivative(k,5)[1] = -x*z;
-              data.derivative(k,6)[1] = (1.-x)*z;
-              data.derivative(k,7)[1] = x*z;
-              data.derivative(k,0)[2] = (x-1)*(1.-y);
-              data.derivative(k,1)[2] = x*(y-1.);
-              data.derivative(k,2)[2] = (x-1.)*y;
-              data.derivative(k,3)[2] = -x*y;
-              data.derivative(k,4)[2] = (1.-x)*(1.-y);
-              data.derivative(k,5)[2] = x*(1.-y);
-              data.derivative(k,6)[2] = (1.-x)*y;
-              data.derivative(k,7)[2] = x*y;
-            }
-          if (data.shape_second_derivatives.size()!=0)
-            {
-              // the following may or may not
-              // work if dim != spacedim
-              Assert (spacedim == 3, ExcNotImplemented());
-
-              Assert(data.shape_second_derivatives.size()==n_shape_functions*n_points,
-                     ExcInternalError());
-              data.second_derivative(k,0)[0][0] = 0;
-              data.second_derivative(k,1)[0][0] = 0;
-              data.second_derivative(k,2)[0][0] = 0;
-              data.second_derivative(k,3)[0][0] = 0;
-              data.second_derivative(k,4)[0][0] = 0;
-              data.second_derivative(k,5)[0][0] = 0;
-              data.second_derivative(k,6)[0][0] = 0;
-              data.second_derivative(k,7)[0][0] = 0;
-              data.second_derivative(k,0)[1][1] = 0;
-              data.second_derivative(k,1)[1][1] = 0;
-              data.second_derivative(k,2)[1][1] = 0;
-              data.second_derivative(k,3)[1][1] = 0;
-              data.second_derivative(k,4)[1][1] = 0;
-              data.second_derivative(k,5)[1][1] = 0;
-              data.second_derivative(k,6)[1][1] = 0;
-              data.second_derivative(k,7)[1][1] = 0;
-              data.second_derivative(k,0)[2][2] = 0;
-              data.second_derivative(k,1)[2][2] = 0;
-              data.second_derivative(k,2)[2][2] = 0;
-              data.second_derivative(k,3)[2][2] = 0;
-              data.second_derivative(k,4)[2][2] = 0;
-              data.second_derivative(k,5)[2][2] = 0;
-              data.second_derivative(k,6)[2][2] = 0;
-              data.second_derivative(k,7)[2][2] = 0;
-
-              data.second_derivative(k,0)[0][1] = (1.-z);
-              data.second_derivative(k,1)[0][1] = -(1.-z);
-              data.second_derivative(k,2)[0][1] = -(1.-z);
-              data.second_derivative(k,3)[0][1] = (1.-z);
-              data.second_derivative(k,4)[0][1] = z;
-              data.second_derivative(k,5)[0][1] = -z;
-              data.second_derivative(k,6)[0][1] = -z;
-              data.second_derivative(k,7)[0][1] = z;
-              data.second_derivative(k,0)[1][0] = (1.-z);
-              data.second_derivative(k,1)[1][0] = -(1.-z);
-              data.second_derivative(k,2)[1][0] = -(1.-z);
-              data.second_derivative(k,3)[1][0] = (1.-z);
-              data.second_derivative(k,4)[1][0] = z;
-              data.second_derivative(k,5)[1][0] = -z;
-              data.second_derivative(k,6)[1][0] = -z;
-              data.second_derivative(k,7)[1][0] = z;
-
-              data.second_derivative(k,0)[0][2] = (1.-y);
-              data.second_derivative(k,1)[0][2] = -(1.-y);
-              data.second_derivative(k,2)[0][2] = y;
-              data.second_derivative(k,3)[0][2] = -y;
-              data.second_derivative(k,4)[0][2] = -(1.-y);
-              data.second_derivative(k,5)[0][2] = (1.-y);
-              data.second_derivative(k,6)[0][2] = -y;
-              data.second_derivative(k,7)[0][2] = y;
-              data.second_derivative(k,0)[2][0] = (1.-y);
-              data.second_derivative(k,1)[2][0] = -(1.-y);
-              data.second_derivative(k,2)[2][0] = y;
-              data.second_derivative(k,3)[2][0] = -y;
-              data.second_derivative(k,4)[2][0] = -(1.-y);
-              data.second_derivative(k,5)[2][0] = (1.-y);
-              data.second_derivative(k,6)[2][0] = -y;
-              data.second_derivative(k,7)[2][0] = y;
-
-              data.second_derivative(k,0)[1][2] = (1.-x);
-              data.second_derivative(k,1)[1][2] = x;
-              data.second_derivative(k,2)[1][2] = -(1.-x);
-              data.second_derivative(k,3)[1][2] = -x;
-              data.second_derivative(k,4)[1][2] = -(1.-x);
-              data.second_derivative(k,5)[1][2] = -x;
-              data.second_derivative(k,6)[1][2] = (1.-x);
-              data.second_derivative(k,7)[1][2] = x;
-              data.second_derivative(k,0)[2][1] = (1.-x);
-              data.second_derivative(k,1)[2][1] = x;
-              data.second_derivative(k,2)[2][1] = -(1.-x);
-              data.second_derivative(k,3)[2][1] = -x;
-              data.second_derivative(k,4)[2][1] = -(1.-x);
-              data.second_derivative(k,5)[2][1] = -x;
-              data.second_derivative(k,6)[2][1] = (1.-x);
-              data.second_derivative(k,7)[2][1] = x;
-            }
-          if (data.shape_third_derivatives.size()!=0)
-            {
-              // if lower order derivative don't work, neither should this
-              Assert (spacedim == 3, ExcNotImplemented());
-
-              Assert(data.shape_third_derivatives.size()==n_shape_functions*n_points,
-                     ExcInternalError());
-
-              for (unsigned int i=0; i<3; ++i)
-                for (unsigned int j=0; j<3; ++j)
-                  for (unsigned int l=0; l<3; ++l)
-                    if ((i==j)||(j==l)||(l==i))
-                      {
-                        for (unsigned int m=0; m<8; ++m)
-                          data.third_derivative(k,m)[i][j][l] = 0;
-                      }
-                    else
-                      {
-                        data.third_derivative(k,0)[i][j][l] = -1.;
-                        data.third_derivative(k,1)[i][j][l] = 1.;
-                        data.third_derivative(k,2)[i][j][l] = 1.;
-                        data.third_derivative(k,3)[i][j][l] = -1.;
-                        data.third_derivative(k,4)[i][j][l] = 1.;
-                        data.third_derivative(k,5)[i][j][l] = -1.;
-                        data.third_derivative(k,6)[i][j][l] = -1.;
-                        data.third_derivative(k,7)[i][j][l] = 1.;
-                      }
-
-            }
-          if (data.shape_fourth_derivatives.size()!=0)
-            {
-              // if lower order derivative don't work, neither should this
-              Assert (spacedim == 3, ExcNotImplemented());
-
-              Assert(data.shape_fourth_derivatives.size()==n_shape_functions*n_points,
-                     ExcInternalError());
-              Tensor<4,3> zero;
-              for (unsigned int i=0; i<8; ++i)
-                data.fourth_derivative(k,i) = zero;
+                             unit_tangentials[nfaces+i].end(), tang2);
+                }
             }
         }
     }
-  }
 }
 
 
+
 namespace
 {
   template <int dim>
@@ -555,7 +840,7 @@ compute_shape_function_values (const std::vector<Point<dim> > &unit_points)
   if ((polynomial_degree == 1)
       &&
       (dim == spacedim))
-    internal::MappingQGeneric::compute_shape_function_values<spacedim> (n_shape_functions,
+    internal::MappingQ1::compute_shape_function_values<spacedim> (n_shape_functions,
         unit_points, *this);
   else
     // otherwise ask an object that describes the polynomial space
@@ -1136,595 +1421,303 @@ namespace
 
         ++newton_iteration;
         if (newton_iteration > newton_iteration_limit)
-          AssertThrow (false,
-                       (typename Mapping<dim,spacedim>::ExcTransformationFailed()));
-        last_f_weighted_norm = (df_inverse * f).norm();
-      }
-    while (last_f_weighted_norm > eps);
-
-    return p_unit;
-  }
-
-
-
-  /**
-   * Implementation of transform_real_to_unit_cell for dim==spacedim-1
-   */
-  template <int dim>
-  Point<dim>
-  do_transform_real_to_unit_cell_internal_codim1
-  (const typename Triangulation<dim,dim+1>::cell_iterator &cell,
-   const Point<dim+1>                                       &p,
-   const Point<dim>                                         &initial_p_unit,
-   typename MappingQGeneric<dim,dim+1>::InternalData       &mdata)
-  {
-    const unsigned int spacedim = dim+1;
-
-    const unsigned int n_shapes=mdata.shape_values.size();
-    (void)n_shapes;
-    Assert(n_shapes!=0, ExcInternalError());
-    Assert(mdata.shape_derivatives.size()==n_shapes, ExcInternalError());
-    Assert(mdata.shape_second_derivatives.size()==n_shapes, ExcInternalError());
-
-    std::vector<Point<spacedim> > &points=mdata.mapping_support_points;
-    Assert(points.size()==n_shapes, ExcInternalError());
-
-    Point<spacedim> p_minus_F;
-
-    Tensor<1,spacedim>  DF[dim];
-    Tensor<1,spacedim>  D2F[dim][dim];
-
-    Point<dim> p_unit = initial_p_unit;
-    Point<dim> f;
-    Tensor<2,dim>  df;
-
-    // Evaluate first and second derivatives
-    mdata.compute_shape_function_values(std::vector<Point<dim> > (1, p_unit));
-
-    for (unsigned int k=0; k<mdata.n_shape_functions; ++k)
-      {
-        const Tensor<1,dim>   &grad_phi_k = mdata.derivative(0,k);
-        const Tensor<2,dim>   &hessian_k  = mdata.second_derivative(0,k);
-        const Point<spacedim> &point_k = points[k];
-
-        for (unsigned int j=0; j<dim; ++j)
-          {
-            DF[j] += grad_phi_k[j] * point_k;
-            for (unsigned int l=0; l<dim; ++l)
-              D2F[j][l] += hessian_k[j][l] * point_k;
-          }
-      }
-
-    p_minus_F = p;
-    p_minus_F -= compute_mapped_location_of_point<dim,spacedim>(mdata);
-
-
-    for (unsigned int j=0; j<dim; ++j)
-      f[j] = DF[j] * p_minus_F;
-
-    for (unsigned int j=0; j<dim; ++j)
-      {
-        f[j] = DF[j] * p_minus_F;
-        for (unsigned int l=0; l<dim; ++l)
-          df[j][l] = -DF[j]*DF[l] + D2F[j][l] * p_minus_F;
-      }
-
-
-    const double eps = 1.e-12*cell->diameter();
-    const unsigned int loop_limit = 10;
-
-    unsigned int loop=0;
-
-    while (f.norm()>eps && loop++<loop_limit)
-      {
-        // Solve  [df(x)]d=f(x)
-        const Tensor<1,dim> d = invert(df) * static_cast<const Tensor<1,dim>&>(f);
-        p_unit -= d;
-
-        for (unsigned int j=0; j<dim; ++j)
-          {
-            DF[j].clear();
-            for (unsigned int l=0; l<dim; ++l)
-              D2F[j][l].clear();
-          }
-
-        mdata.compute_shape_function_values(std::vector<Point<dim> > (1, p_unit));
-
-        for (unsigned int k=0; k<mdata.n_shape_functions; ++k)
-          {
-            const Tensor<1,dim>   &grad_phi_k = mdata.derivative(0,k);
-            const Tensor<2,dim>   &hessian_k  = mdata.second_derivative(0,k);
-            const Point<spacedim> &point_k = points[k];
-
-            for (unsigned int j=0; j<dim; ++j)
-              {
-                DF[j] += grad_phi_k[j] * point_k;
-                for (unsigned int l=0; l<dim; ++l)
-                  D2F[j][l] += hessian_k[j][l] * point_k;
-              }
-          }
-
-        //TODO: implement a line search here in much the same way as for
-        // the corresponding function above that does so for dim==spacedim
-        p_minus_F = p;
-        p_minus_F -= compute_mapped_location_of_point<dim,spacedim>(mdata);
-
-        for (unsigned int j=0; j<dim; ++j)
-          {
-            f[j] = DF[j] * p_minus_F;
-            for (unsigned int l=0; l<dim; ++l)
-              df[j][l] = -DF[j]*DF[l] + D2F[j][l] * p_minus_F;
-          }
-
-      }
-
-
-    // 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 p_unit;
-  }
-
-
-}
-
-
-
-// visual studio freaks out when trying to determine if
-// do_transform_real_to_unit_cell_internal with dim=3 and spacedim=4 is a good
-// candidate. So instead of letting the compiler pick the correct overload, we
-// use template specialization to make sure we pick up the right function to
-// call:
-
-template<int dim, int spacedim>
-Point<dim>
-MappingQGeneric<dim,spacedim>::
-transform_real_to_unit_cell_internal
-(const typename Triangulation<dim,spacedim>::cell_iterator &,
- const Point<spacedim> &,
- const Point<dim> &) const
-{
-  // default implementation (should never be called)
-  Assert(false, ExcInternalError());
-  return Point<dim>();
-}
-
-template<>
-Point<1>
-MappingQGeneric<1,1>::
-transform_real_to_unit_cell_internal
-(const Triangulation<1,1>::cell_iterator &cell,
- const Point<1>                            &p,
- const Point<1>                                 &initial_p_unit) const
-{
-  const int dim = 1;
-  const int spacedim = 1;
-
-  const Quadrature<dim> point_quadrature(initial_p_unit);
-
-  UpdateFlags update_flags = update_quadrature_points | update_jacobians;
-  if (spacedim>dim)
-    update_flags |= update_jacobian_grads;
-  std_cxx11::unique_ptr<InternalData> mdata (get_data(update_flags,
-                                                      point_quadrature));
-
-  mdata->mapping_support_points = this->compute_mapping_support_points (cell);
-
-  // dispatch to the various specializations for spacedim=dim,
-  // spacedim=dim+1, etc
-  return do_transform_real_to_unit_cell_internal<1>(cell, p, initial_p_unit, *mdata);
-}
-
-template<>
-Point<2>
-MappingQGeneric<2, 2>::
-transform_real_to_unit_cell_internal
-(const Triangulation<2, 2>::cell_iterator &cell,
- const Point<2>                            &p,
- const Point<2>                                 &initial_p_unit) const
-{
-  const int dim = 2;
-  const int spacedim = 2;
-
-  const Quadrature<dim> point_quadrature(initial_p_unit);
-
-  UpdateFlags update_flags = update_quadrature_points | update_jacobians;
-  if (spacedim>dim)
-    update_flags |= update_jacobian_grads;
-  std_cxx11::unique_ptr<InternalData> mdata (get_data(update_flags,
-                                                      point_quadrature));
-
-  mdata->mapping_support_points = this->compute_mapping_support_points (cell);
-
-  // dispatch to the various specializations for spacedim=dim,
-  // spacedim=dim+1, etc
-  return do_transform_real_to_unit_cell_internal<2>(cell, p, initial_p_unit, *mdata);
-}
-
-template<>
-Point<3>
-MappingQGeneric<3, 3>::
-transform_real_to_unit_cell_internal
-(const Triangulation<3, 3>::cell_iterator &cell,
- const Point<3>                            &p,
- const Point<3>                                 &initial_p_unit) const
-{
-  const int dim = 3;
-  const int spacedim = 3;
-
-  const Quadrature<dim> point_quadrature(initial_p_unit);
-
-  UpdateFlags update_flags = update_quadrature_points | update_jacobians;
-  if (spacedim>dim)
-    update_flags |= update_jacobian_grads;
-  std_cxx11::unique_ptr<InternalData> mdata (get_data(update_flags,
-                                                      point_quadrature));
-
-  mdata->mapping_support_points = this->compute_mapping_support_points (cell);
-
-  // dispatch to the various specializations for spacedim=dim,
-  // spacedim=dim+1, etc
-  return do_transform_real_to_unit_cell_internal<3>(cell, p, initial_p_unit, *mdata);
-}
-
-template<>
-Point<1>
-MappingQGeneric<1, 2>::
-transform_real_to_unit_cell_internal
-(const Triangulation<1, 2>::cell_iterator &cell,
- const Point<2>                            &p,
- const Point<1>                                 &initial_p_unit) const
-{
-  const int dim = 1;
-  const int spacedim = 2;
-
-  const Quadrature<dim> point_quadrature(initial_p_unit);
-
-  UpdateFlags update_flags = update_quadrature_points | update_jacobians;
-  if (spacedim>dim)
-    update_flags |= update_jacobian_grads;
-  std_cxx11::unique_ptr<InternalData> mdata (get_data(update_flags,
-                                                      point_quadrature));
+          AssertThrow (false,
+                       (typename Mapping<dim,spacedim>::ExcTransformationFailed()));
+        last_f_weighted_norm = (df_inverse * f).norm();
+      }
+    while (last_f_weighted_norm > eps);
 
-  mdata->mapping_support_points = this->compute_mapping_support_points (cell);
+    return p_unit;
+  }
 
-  // dispatch to the various specializations for spacedim=dim,
-  // spacedim=dim+1, etc
-  return do_transform_real_to_unit_cell_internal_codim1<1>(cell, p, initial_p_unit, *mdata);
-}
 
-template<>
-Point<2>
-MappingQGeneric<2, 3>::
-transform_real_to_unit_cell_internal
-(const Triangulation<2, 3>::cell_iterator &cell,
- const Point<3>                            &p,
- const Point<2>                                 &initial_p_unit) const
-{
-  const int dim = 2;
-  const int spacedim = 3;
 
-  const Quadrature<dim> point_quadrature(initial_p_unit);
+  /**
+   * Implementation of transform_real_to_unit_cell for dim==spacedim-1
+   */
+  template <int dim>
+  Point<dim>
+  do_transform_real_to_unit_cell_internal_codim1
+  (const typename Triangulation<dim,dim+1>::cell_iterator &cell,
+   const Point<dim+1>                                       &p,
+   const Point<dim>                                         &initial_p_unit,
+   typename MappingQGeneric<dim,dim+1>::InternalData       &mdata)
+  {
+    const unsigned int spacedim = dim+1;
 
-  UpdateFlags update_flags = update_quadrature_points | update_jacobians;
-  if (spacedim>dim)
-    update_flags |= update_jacobian_grads;
-  std_cxx11::unique_ptr<InternalData> mdata (get_data(update_flags,
-                                                      point_quadrature));
+    const unsigned int n_shapes=mdata.shape_values.size();
+    (void)n_shapes;
+    Assert(n_shapes!=0, ExcInternalError());
+    Assert(mdata.shape_derivatives.size()==n_shapes, ExcInternalError());
+    Assert(mdata.shape_second_derivatives.size()==n_shapes, ExcInternalError());
 
-  mdata->mapping_support_points = this->compute_mapping_support_points (cell);
+    std::vector<Point<spacedim> > &points=mdata.mapping_support_points;
+    Assert(points.size()==n_shapes, ExcInternalError());
 
-  // dispatch to the various specializations for spacedim=dim,
-  // spacedim=dim+1, etc
-  return do_transform_real_to_unit_cell_internal_codim1<2>(cell, p, initial_p_unit, *mdata);
-}
+    Point<spacedim> p_minus_F;
 
-template<>
-Point<1>
-MappingQGeneric<1, 3>::
-transform_real_to_unit_cell_internal
-(const Triangulation<1, 3>::cell_iterator &,
- const Point<3> &,
- const Point<1> &) const
-{
-  Assert (false, ExcNotImplemented());
-  return Point<1>();
-}
+    Tensor<1,spacedim>  DF[dim];
+    Tensor<1,spacedim>  D2F[dim][dim];
 
+    Point<dim> p_unit = initial_p_unit;
+    Point<dim> f;
+    Tensor<2,dim>  df;
 
-namespace internal
-{
-  namespace MappingQ1
-  {
-    namespace
-    {
+    // Evaluate first and second derivatives
+    mdata.compute_shape_function_values(std::vector<Point<dim> > (1, p_unit));
 
-      // 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)
+    for (unsigned int k=0; k<mdata.n_shape_functions; ++k)
       {
-        Assert(spacedim == 1, ExcInternalError());
-        return Point<1>((p[0] - vertices[0](0))/(vertices[1](0) - vertices[0](0)));
+        const Tensor<1,dim>   &grad_phi_k = mdata.derivative(0,k);
+        const Tensor<2,dim>   &hessian_k  = mdata.second_derivative(0,k);
+        const Point<spacedim> &point_k = points[k];
+
+        for (unsigned int j=0; j<dim; ++j)
+          {
+            DF[j] += grad_phi_k[j] * point_k;
+            for (unsigned int l=0; l<dim; ++l)
+              D2F[j][l] += hessian_k[j][l] * point_k;
+          }
       }
 
+    p_minus_F = p;
+    p_minus_F -= compute_mapped_location_of_point<dim,spacedim>(mdata);
 
 
-      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)
+    for (unsigned int j=0; j<dim; ++j)
+      f[j] = DF[j] * p_minus_F;
+
+    for (unsigned int j=0; j<dim; ++j)
       {
-        Assert(spacedim == 2, ExcInternalError());
-        const double x = p(0);
-        const double y = p(1);
+        f[j] = DF[j] * p_minus_F;
+        for (unsigned int l=0; l<dim; ++l)
+          df[j][l] = -DF[j]*DF[l] + D2F[j][l] * p_minus_F;
+      }
 
-        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 eps = 1.e-12*cell->diameter();
+    const unsigned int loop_limit = 10;
 
-        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;
+    unsigned int loop=0;
 
-        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()));
-          }
+    while (f.norm()>eps && loop++<loop_limit)
+      {
+        // Solve  [df(x)]d=f(x)
+        const Tensor<1,dim> d = invert(df) * static_cast<const Tensor<1,dim>&>(f);
+        p_unit -= d;
 
-        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
+        for (unsigned int j=0; j<dim; ++j)
           {
-            eta1 = 2*c / (-b - std::sqrt(discriminant));
-            eta2 = 2*c / (-b + std::sqrt(discriminant));
+            DF[j].clear();
+            for (unsigned int l=0; l<dim; ++l)
+              D2F[j][l].clear();
           }
-        // 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)));
+        mdata.compute_shape_function_values(std::vector<Point<dim> > (1, p_unit));
 
-        if (std::abs(xi_denominator0) > 1e-10*max_x)
-          {
-            const double xi = (x + subexpr0)/xi_denominator0;
-            return Point<2>(xi, eta);
-          }
-        else
+        for (unsigned int k=0; k<mdata.n_shape_functions; ++k)
           {
-            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
+            const Tensor<1,dim>   &grad_phi_k = mdata.derivative(0,k);
+            const Tensor<2,dim>   &hessian_k  = mdata.second_derivative(0,k);
+            const Point<spacedim> &point_k = points[k];
+
+            for (unsigned int j=0; j<dim; ++j)
               {
-                AssertThrow (false,
-                             (typename Mapping<spacedim,spacedim>::ExcTransformationFailed()));
+                DF[j] += grad_phi_k[j] * point_k;
+                for (unsigned int l=0; l<dim; ++l)
+                  D2F[j][l] += hessian_k[j][l] * point_k;
               }
           }
-        // 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());
-      }
 
+        //TODO: implement a line search here in much the same way as for
+        // the corresponding function above that does so for dim==spacedim
+        p_minus_F = p;
+        p_minus_F -= compute_mapped_location_of_point<dim,spacedim>(mdata);
 
+        for (unsigned int j=0; j<dim; ++j)
+          {
+            f[j] = DF[j] * p_minus_F;
+            for (unsigned int l=0; l<dim; ++l)
+              df[j][l] = -DF[j]*DF[l] + D2F[j][l] * p_minus_F;
+          }
 
-      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>();
       }
 
 
+    // 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 p_unit;
+  }
+
+
+}
+
+
+
+// visual studio freaks out when trying to determine if
+// do_transform_real_to_unit_cell_internal with dim=3 and spacedim=4 is a good
+// candidate. So instead of letting the compiler pick the correct overload, we
+// use template specialization to make sure we pick up the right function to
+// call:
+
+template<int dim, int spacedim>
+Point<dim>
+MappingQGeneric<dim,spacedim>::
+transform_real_to_unit_cell_internal
+(const typename Triangulation<dim,spacedim>::cell_iterator &,
+ const Point<spacedim> &,
+ const Point<dim> &) const
+{
+  // default implementation (should never be called)
+  Assert(false, ExcInternalError());
+  return Point<dim>();
+}
+
+template<>
+Point<1>
+MappingQGeneric<1,1>::
+transform_real_to_unit_cell_internal
+(const Triangulation<1,1>::cell_iterator &cell,
+ const Point<1>                            &p,
+ const Point<1>                                 &initial_p_unit) const
+{
+  const int dim = 1;
+  const int spacedim = 1;
+
+  const Quadrature<dim> point_quadrature(initial_p_unit);
 
-      /**
-       * 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];
-      };
+  UpdateFlags update_flags = update_quadrature_points | update_jacobians;
+  if (spacedim>dim)
+    update_flags |= update_jacobian_grads;
+  std_cxx11::unique_ptr<InternalData> mdata (get_data(update_flags,
+                                                      point_quadrature));
 
+  mdata->mapping_support_points = this->compute_mapping_support_points (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}
-      };
+  // dispatch to the various specializations for spacedim=dim,
+  // spacedim=dim+1, etc
+  return do_transform_real_to_unit_cell_internal<1>(cell, p, initial_p_unit, *mdata);
+}
 
-      template <>
-      const double
-      TransformR2UInitialGuess<1>::
-      Kb[GeometryInfo<1>::vertices_per_cell] = {1.000000, 0.000000};
+template<>
+Point<2>
+MappingQGeneric<2, 2>::
+transform_real_to_unit_cell_internal
+(const Triangulation<2, 2>::cell_iterator &cell,
+ const Point<2>                            &p,
+ const Point<2>                                 &initial_p_unit) const
+{
+  const int dim = 2;
+  const int spacedim = 2;
 
+  const Quadrature<dim> point_quadrature(initial_p_unit);
 
-      /*
-        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}
-      };
+  UpdateFlags update_flags = update_quadrature_points | update_jacobians;
+  if (spacedim>dim)
+    update_flags |= update_jacobian_grads;
+  std_cxx11::unique_ptr<InternalData> mdata (get_data(update_flags,
+                                                      point_quadrature));
 
-      /*
-        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 };
+  mdata->mapping_support_points = this->compute_mapping_support_points (cell);
 
+  // dispatch to the various specializations for spacedim=dim,
+  // spacedim=dim+1, etc
+  return do_transform_real_to_unit_cell_internal<2>(cell, p, initial_p_unit, *mdata);
+}
 
-      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<>
+Point<3>
+MappingQGeneric<3, 3>::
+transform_real_to_unit_cell_internal
+(const Triangulation<3, 3>::cell_iterator &cell,
+ const Point<3>                            &p,
+ const Point<3>                                 &initial_p_unit) const
+{
+  const int dim = 3;
+  const int spacedim = 3;
 
-      };
+  const Quadrature<dim> point_quadrature(initial_p_unit);
 
+  UpdateFlags update_flags = update_quadrature_points | update_jacobians;
+  if (spacedim>dim)
+    update_flags |= update_jacobian_grads;
+  std_cxx11::unique_ptr<InternalData> mdata (get_data(update_flags,
+                                                      point_quadrature));
 
-      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};
+  mdata->mapping_support_points = this->compute_mapping_support_points (cell);
 
-      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;
+  // dispatch to the various specializations for spacedim=dim,
+  // spacedim=dim+1, etc
+  return do_transform_real_to_unit_cell_internal<3>(cell, p, initial_p_unit, *mdata);
+}
 
-        dealii::FullMatrix<double>  KA(GeometryInfo<dim>::vertices_per_cell, dim);
-        dealii::Vector <double>  Kb(GeometryInfo<dim>::vertices_per_cell);
+template<>
+Point<1>
+MappingQGeneric<1, 2>::
+transform_real_to_unit_cell_internal
+(const Triangulation<1, 2>::cell_iterator &cell,
+ const Point<2>                            &p,
+ const Point<1>                                 &initial_p_unit) const
+{
+  const int dim = 1;
+  const int spacedim = 2;
 
-        KA.fill( (double *)(TransformR2UInitialGuess<dim>::KA) );
-        for (unsigned int i=0; i<GeometryInfo<dim>::vertices_per_cell; ++i)
-          Kb(i) = TransformR2UInitialGuess<dim>::Kb[i];
+  const Quadrature<dim> point_quadrature(initial_p_unit);
 
-        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];
+  UpdateFlags update_flags = update_quadrature_points | update_jacobians;
+  if (spacedim>dim)
+    update_flags |= update_jacobian_grads;
+  std_cxx11::unique_ptr<InternalData> mdata (get_data(update_flags,
+                                                      point_quadrature));
 
-        FullMatrix<double> A(spacedim,dim);
-        Y.mmult(A,KA); // A = Y*KA
-        dealii::Vector<double> b(spacedim);
-        Y.vmult(b,Kb); // b = Y*Kb
+  mdata->mapping_support_points = this->compute_mapping_support_points (cell);
 
-        for (unsigned int i=0; i<spacedim; ++i)
-          b(i) -= p[i];
-        b*=-1;
+  // dispatch to the various specializations for spacedim=dim,
+  // spacedim=dim+1, etc
+  return do_transform_real_to_unit_cell_internal_codim1<1>(cell, p, initial_p_unit, *mdata);
+}
 
-        dealii::Vector<double> dest(dim);
+template<>
+Point<2>
+MappingQGeneric<2, 3>::
+transform_real_to_unit_cell_internal
+(const Triangulation<2, 3>::cell_iterator &cell,
+ const Point<3>                            &p,
+ const Point<2>                                 &initial_p_unit) const
+{
+  const int dim = 2;
+  const int spacedim = 3;
 
-        FullMatrix<double> A_1(dim,spacedim);
-        if (dim<spacedim)
-          A_1.left_invert(A);
-        else
-          A_1.invert(A);
+  const Quadrature<dim> point_quadrature(initial_p_unit);
 
-        A_1.vmult(dest,b); //A^{-1}*b
+  UpdateFlags update_flags = update_quadrature_points | update_jacobians;
+  if (spacedim>dim)
+    update_flags |= update_jacobian_grads;
+  std_cxx11::unique_ptr<InternalData> mdata (get_data(update_flags,
+                                                      point_quadrature));
 
-        for (unsigned int i=0; i<dim; ++i)
-          p_unit[i]=dest(i);
+  mdata->mapping_support_points = this->compute_mapping_support_points (cell);
 
-        return p_unit;
-      }
-    }
-  }
+  // dispatch to the various specializations for spacedim=dim,
+  // spacedim=dim+1, etc
+  return do_transform_real_to_unit_cell_internal_codim1<2>(cell, p, initial_p_unit, *mdata);
 }
 
+template<>
+Point<1>
+MappingQGeneric<1, 3>::
+transform_real_to_unit_cell_internal
+(const Triangulation<1, 3>::cell_iterator &,
+ const Point<3> &,
+ const Point<1> &) const
+{
+  Assert (false, ExcNotImplemented());
+  return Point<1>();
+}
 
 
 

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.