From ab3335d8c7a1ccf618eee7f0a5df3d5ea8a3a70c Mon Sep 17 00:00:00 2001 From: David Wells Date: Tue, 20 Oct 2015 11:15:08 -0400 Subject: [PATCH] Remove the internal::MappingQGeneric namespace. 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 | 1723 ++++++++++++++++---------------- 1 file changed, 858 insertions(+), 865 deletions(-) diff --git a/source/fe/mapping_q_generic.cc b/source/fe/mapping_q_generic.cc index 91db94cfc0..7095eba7e9 100644 --- a/source/fe/mapping_q_generic.cc +++ b/source/fe/mapping_q_generic.cc @@ -40,6 +40,624 @@ 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 + Point<1> + transform_real_to_unit_cell + (const std_cxx11::array, GeometryInfo<1>::vertices_per_cell> &vertices, + const Point &p) + { + Assert(spacedim == 1, ExcInternalError()); + return Point<1>((p[0] - vertices[0](0))/(vertices[1](0) - vertices[0](0))); + } + + + + template + Point<2> + transform_real_to_unit_cell + (const std_cxx11::array, GeometryInfo<2>::vertices_per_cell> &vertices, + const Point &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::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::ExcTransformationFailed())); + } + } + // bogus return to placate compiler. It should not be possible to get + // here. + Assert(false, ExcInternalError()); + return Point<2>(std::numeric_limits::quiet_NaN(), + std::numeric_limits::quiet_NaN()); + } + + + + template + Point<3> + transform_real_to_unit_cell + (const std_cxx11::array, GeometryInfo<3>::vertices_per_cell> &/*vertices*/, + const Point &/*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: + *
    + *
  • 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. + *
  • Solve the equation A x_hat + b = p for x_hat + *
  • This x_hat is the initial solution used for the Newton Method. + *
+ * + * @note if dim + struct TransformR2UInitialGuess + { + static const double KA[GeometryInfo::vertices_per_cell][dim]; + static const double Kb[GeometryInfo::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 + Point + transform_real_to_unit_cell_initial_guess (const std::vector > &vertex, + const Point &p) + { + Point p_unit; + + dealii::FullMatrix KA(GeometryInfo::vertices_per_cell, dim); + dealii::Vector Kb(GeometryInfo::vertices_per_cell); + + KA.fill( (double *)(TransformR2UInitialGuess::KA) ); + for (unsigned int i=0; i::vertices_per_cell; ++i) + Kb(i) = TransformR2UInitialGuess::Kb[i]; + + FullMatrix Y(spacedim, GeometryInfo::vertices_per_cell); + for (unsigned int v=0; v::vertices_per_cell; v++) + for (unsigned int i=0; i A(spacedim,dim); + Y.mmult(A,KA); // A = Y*KA + dealii::Vector b(spacedim); + Y.vmult(b,Kb); // b = Y*Kb + + for (unsigned int i=0; i dest(dim); + + FullMatrix A_1(dim,spacedim); + if (dim + void + compute_shape_function_values (const unsigned int n_shape_functions, + const std::vector > &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 + void + compute_shape_function_values (const unsigned int n_shape_functions, + const std::vector > &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 + void + compute_shape_function_values (const unsigned int n_shape_functions, + const std::vector > &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 MappingQGeneric::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 - void - compute_shape_function_values (const unsigned int n_shape_functions, - const std::vector > &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 - void - compute_shape_function_values (const unsigned int n_shape_functions, - const std::vector > &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 - void - compute_shape_function_values (const unsigned int n_shape_functions, - const std::vector > &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 @@ -555,7 +840,7 @@ compute_shape_function_values (const std::vector > &unit_points) if ((polynomial_degree == 1) && (dim == spacedim)) - internal::MappingQGeneric::compute_shape_function_values (n_shape_functions, + internal::MappingQ1::compute_shape_function_values (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::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 - Point - do_transform_real_to_unit_cell_internal_codim1 - (const typename Triangulation::cell_iterator &cell, - const Point &p, - const Point &initial_p_unit, - typename MappingQGeneric::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 > &points=mdata.mapping_support_points; - Assert(points.size()==n_shapes, ExcInternalError()); - - Point p_minus_F; - - Tensor<1,spacedim> DF[dim]; - Tensor<1,spacedim> D2F[dim][dim]; - - Point p_unit = initial_p_unit; - Point f; - Tensor<2,dim> df; - - // Evaluate first and second derivatives - mdata.compute_shape_function_values(std::vector > (1, p_unit)); - - for (unsigned int k=0; k &grad_phi_k = mdata.derivative(0,k); - const Tensor<2,dim> &hessian_k = mdata.second_derivative(0,k); - const Point &point_k = points[k]; - - for (unsigned int j=0; j(mdata); - - - for (unsigned int j=0; jdiameter(); - const unsigned int loop_limit = 10; - - unsigned int loop=0; - - while (f.norm()>eps && loop++ d = invert(df) * static_cast&>(f); - p_unit -= d; - - for (unsigned int j=0; j > (1, p_unit)); - - for (unsigned int k=0; k &grad_phi_k = mdata.derivative(0,k); - const Tensor<2,dim> &hessian_k = mdata.second_derivative(0,k); - const Point &point_k = points[k]; - - for (unsigned int j=0; j(mdata); - - for (unsigned int j=0; j::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 -Point -MappingQGeneric:: -transform_real_to_unit_cell_internal -(const typename Triangulation::cell_iterator &, - const Point &, - const Point &) const -{ - // default implementation (should never be called) - Assert(false, ExcInternalError()); - return Point(); -} - -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 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 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 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 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 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 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 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 mdata (get_data(update_flags, - point_quadrature)); + AssertThrow (false, + (typename Mapping::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 point_quadrature(initial_p_unit); + /** + * Implementation of transform_real_to_unit_cell for dim==spacedim-1 + */ + template + Point + do_transform_real_to_unit_cell_internal_codim1 + (const typename Triangulation::cell_iterator &cell, + const Point &p, + const Point &initial_p_unit, + typename MappingQGeneric::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 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 > &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 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 p_unit = initial_p_unit; + Point f; + Tensor<2,dim> df; -namespace internal -{ - namespace MappingQ1 - { - namespace - { + // Evaluate first and second derivatives + mdata.compute_shape_function_values(std::vector > (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 - Point<1> - transform_real_to_unit_cell - (const std_cxx11::array, GeometryInfo<1>::vertices_per_cell> &vertices, - const Point &p) + for (unsigned int k=0; k((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 &point_k = points[k]; + + for (unsigned int j=0; j(mdata); - template - Point<2> - transform_real_to_unit_cell - (const std_cxx11::array, GeometryInfo<2>::vertices_per_cell> &vertices, - const Point &p) + for (unsigned int j=0; jdiameter(); + 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::ExcTransformationFailed())); - } + while (f.norm()>eps && loop++ d = invert(df) * static_cast&>(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 > (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 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 &point_k = points[k]; + + for (unsigned int j=0; j::ExcTransformationFailed())); + DF[j] += grad_phi_k[j] * point_k; + for (unsigned int l=0; l(std::numeric_limits::quiet_NaN(), - std::numeric_limits::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(mdata); + for (unsigned int j=0; j - Point<3> - transform_real_to_unit_cell - (const std_cxx11::array, GeometryInfo<3>::vertices_per_cell> &/*vertices*/, - const Point &/*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::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 +Point +MappingQGeneric:: +transform_real_to_unit_cell_internal +(const typename Triangulation::cell_iterator &, + const Point &, + const Point &) const +{ + // default implementation (should never be called) + Assert(false, ExcInternalError()); + return Point(); +} + +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 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: - *
    - *
  • 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. - *
  • Solve the equation A x_hat + b = p for x_hat - *
  • This x_hat is the initial solution used for the Newton Method. - *
- * - * @note if dim - struct TransformR2UInitialGuess - { - static const double KA[GeometryInfo::vertices_per_cell][dim]; - static const double Kb[GeometryInfo::vertices_per_cell]; - }; + UpdateFlags update_flags = update_quadrature_points | update_jacobians; + if (spacedim>dim) + update_flags |= update_jacobian_grads; + std_cxx11::unique_ptr 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 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 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 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 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 - Point - transform_real_to_unit_cell_initial_guess (const std::vector > &vertex, - const Point &p) - { - Point 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 KA(GeometryInfo::vertices_per_cell, dim); - dealii::Vector Kb(GeometryInfo::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::KA) ); - for (unsigned int i=0; i::vertices_per_cell; ++i) - Kb(i) = TransformR2UInitialGuess::Kb[i]; + const Quadrature point_quadrature(initial_p_unit); - FullMatrix Y(spacedim, GeometryInfo::vertices_per_cell); - for (unsigned int v=0; v::vertices_per_cell; v++) - for (unsigned int i=0; idim) + update_flags |= update_jacobian_grads; + std_cxx11::unique_ptr mdata (get_data(update_flags, + point_quadrature)); - FullMatrix A(spacedim,dim); - Y.mmult(A,KA); // A = Y*KA - dealii::Vector 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(cell, p, initial_p_unit, *mdata); +} - dealii::Vector 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 A_1(dim,spacedim); - if (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 mdata (get_data(update_flags, + point_quadrature)); - for (unsigned int i=0; imapping_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>(); +} -- 2.39.5