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)
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>
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
++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>();
+}