From de2f4e95cf1bd7c9ba5cf09fc4c8caf8fd90101f Mon Sep 17 00:00:00 2001 From: Martin Kronbichler Date: Wed, 21 Oct 2020 11:30:16 +0200 Subject: [PATCH] Restructure MappingQGeneric: Use internal header file --- include/deal.II/fe/mapping_q_internal.h | 2202 +++++++++++++++++++++++ source/fe/mapping_q_generic.cc | 2199 +--------------------- 2 files changed, 2204 insertions(+), 2197 deletions(-) create mode 100644 include/deal.II/fe/mapping_q_internal.h diff --git a/include/deal.II/fe/mapping_q_internal.h b/include/deal.II/fe/mapping_q_internal.h new file mode 100644 index 0000000000..e0971d8258 --- /dev/null +++ b/include/deal.II/fe/mapping_q_internal.h @@ -0,0 +1,2202 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2020 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- + +#ifndef dealii_mapping_q_internal_h +#define dealii_mapping_q_internal_h + +#include + +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include + +#include +#include +#include +#include + +#include + + +DEAL_II_NAMESPACE_OPEN + +namespace internal +{ + /** + * Internal namespace to implement methods specific to MappingQ1, in + * particular an explicit formula for the transformation from the real to + * the unit cell in 2D. + */ + namespace MappingQ1 + { + // These are left as templates on the spatial dimension (even though dim + // == spacedim must be true for them to make sense) because templates are + // expanded before the compiler eliminates code due to the 'if (dim == + // spacedim)' statement (see the body of the general + // transform_real_to_unit_cell). + template + inline Point<1> + transform_real_to_unit_cell( + const std::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 + inline Point<2> + transform_real_to_unit_cell( + const std::array, GeometryInfo<2>::vertices_per_cell> + & vertices, + const Point &p) + { + Assert(spacedim == 2, ExcInternalError()); + + // For accuracy reasons, we do all arithmetic in extended precision + // (long double). This has a noticeable effect on the hit rate for + // borderline cases and thus makes the algorithm more robust. + const long double x = p(0); + const long double y = p(1); + + const long double x0 = vertices[0](0); + const long double x1 = vertices[1](0); + const long double x2 = vertices[2](0); + const long double x3 = vertices[3](0); + + const long double y0 = vertices[0](1); + const long double y1 = vertices[1](1); + const long double y2 = vertices[2](1); + const long double y3 = vertices[3](1); + + const long double a = (x1 - x3) * (y0 - y2) - (x0 - x2) * (y1 - y3); + const long double b = -(x0 - x1 - x2 + x3) * y + (x - 2 * x1 + x3) * y0 - + (x - 2 * x0 + x2) * y1 - (x - x1) * y2 + + (x - x0) * y3; + const long double c = (x0 - x1) * y - (x - x1) * y0 + (x - x0) * y1; + + const long 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) + AssertThrow( + discriminant > 0.0, + (typename Mapping::ExcTransformationFailed())); + + long double eta1; + long double eta2; + const long double sqrt_discriminant = std::sqrt(discriminant); + // special case #1: if a is near-zero to make the discriminant exactly + // equal b, then use the linear formula + if (b != 0.0 && std::abs(b) == sqrt_discriminant) + { + eta1 = -c / b; + eta2 = -c / b; + } + // special case #2: a is zero for parallelograms and very small for + // near-parallelograms: + else if (std::abs(a) < 1e-8 * std::abs(b)) + { + // if both a and c are very small then the root should be near + // zero: this first case will capture that + eta1 = 2 * c / (-b - sqrt_discriminant); + eta2 = 2 * c / (-b + sqrt_discriminant); + } + // finally, use the plain version: + else + { + eta1 = (-b - sqrt_discriminant) / (2 * a); + eta2 = (-b + sqrt_discriminant) / (2 * a); + } + // pick the one closer to the center of the cell. + const long 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 long double subexpr0 = -eta * x2 + x0 * (eta - 1); + const long double xi_denominator0 = eta * x3 - x1 * (eta - 1) + subexpr0; + const long 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 {xi, static_cast(eta)}; + } + else + { + const long double max_y = + std::max(std::max(std::abs(y0), std::abs(y1)), + std::max(std::abs(y2), std::abs(y3))); + const long double subexpr1 = -eta * y2 + y0 * (eta - 1); + const long 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 {xi, static_cast(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 {std::numeric_limits::quiet_NaN(), + std::numeric_limits::quiet_NaN()}; + } + + + + template + inline Point<3> + transform_real_to_unit_cell( + const std::array, GeometryInfo<3>::vertices_per_cell> + & /*vertices*/, + const Point & /*p*/) + { + // It should not be possible to get here + Assert(false, ExcInternalError()); + return {std::numeric_limits::quiet_NaN(), + std::numeric_limits::quiet_NaN(), + std::numeric_limits::quiet_NaN()}; + } + } // namespace MappingQ1 + + + + /** + * Internal namespace to implement methods of MappingQGeneric, such as the + * evaluation of the mapping and the transformation between real and unit + * cell. + */ + namespace MappingQGenericImplementation + { + /** + * This function is needed by the constructor of + * MappingQ for dim= 2 and 3. + * + * For the definition of the @p support_point_weights_on_quad please + * refer to the description of TransfiniteInterpolationManifold. + */ + inline dealii::Table<2, double> + compute_support_point_weights_on_quad(const unsigned int polynomial_degree) + { + dealii::Table<2, double> loqvs; + + // we are asked to compute weights for interior support points, but + // there are no interior points if degree==1 + if (polynomial_degree == 1) + return loqvs; + + const unsigned int M = polynomial_degree - 1; + const unsigned int n_inner_2d = M * M; + const unsigned int n_outer_2d = 4 + 4 * M; + + // set the weights of transfinite interpolation + loqvs.reinit(n_inner_2d, n_outer_2d); + QGaussLobatto<2> gl(polynomial_degree + 1); + for (unsigned int i = 0; i < M; ++i) + for (unsigned int j = 0; j < M; ++j) + { + const Point<2> p = + gl.point((i + 1) * (polynomial_degree + 1) + (j + 1)); + const unsigned int index_table = i * M + j; + for (unsigned int v = 0; v < 4; ++v) + loqvs(index_table, v) = + -GeometryInfo<2>::d_linear_shape_function(p, v); + loqvs(index_table, 4 + i) = 1. - p[0]; + loqvs(index_table, 4 + i + M) = p[0]; + loqvs(index_table, 4 + j + 2 * M) = 1. - p[1]; + loqvs(index_table, 4 + j + 3 * M) = p[1]; + } + + // the sum of weights of the points at the outer rim should be one. + // check this + for (unsigned int unit_point = 0; unit_point < n_inner_2d; ++unit_point) + Assert(std::fabs(std::accumulate(loqvs[unit_point].begin(), + loqvs[unit_point].end(), + 0.) - + 1) < 1e-13 * polynomial_degree, + ExcInternalError()); + + return loqvs; + } + + + + /** + * This function is needed by the constructor of MappingQ<3>. + * + * For the definition of the @p support_point_weights_on_quad please + * refer to the description of TransfiniteInterpolationManifold. + */ + inline dealii::Table<2, double> + compute_support_point_weights_on_hex(const unsigned int polynomial_degree) + { + dealii::Table<2, double> lohvs; + + // we are asked to compute weights for interior support points, but + // there are no interior points if degree==1 + if (polynomial_degree == 1) + return lohvs; + + const unsigned int M = polynomial_degree - 1; + + const unsigned int n_inner = Utilities::fixed_power<3>(M); + const unsigned int n_outer = 8 + 12 * M + 6 * M * M; + + // set the weights of transfinite interpolation + lohvs.reinit(n_inner, n_outer); + QGaussLobatto<3> gl(polynomial_degree + 1); + for (unsigned int i = 0; i < M; ++i) + for (unsigned int j = 0; j < M; ++j) + for (unsigned int k = 0; k < M; ++k) + { + const Point<3> p = gl.point((i + 1) * (M + 2) * (M + 2) + + (j + 1) * (M + 2) + (k + 1)); + const unsigned int index_table = i * M * M + j * M + k; + + // vertices + for (unsigned int v = 0; v < 8; ++v) + lohvs(index_table, v) = + GeometryInfo<3>::d_linear_shape_function(p, v); + + // lines + { + constexpr std::array line_coordinates_y( + {{0, 1, 4, 5}}); + const Point<2> py(p[0], p[2]); + for (unsigned int l = 0; l < 4; ++l) + lohvs(index_table, 8 + line_coordinates_y[l] * M + j) = + -GeometryInfo<2>::d_linear_shape_function(py, l); + } + + { + constexpr std::array line_coordinates_x( + {{2, 3, 6, 7}}); + const Point<2> px(p[1], p[2]); + for (unsigned int l = 0; l < 4; ++l) + lohvs(index_table, 8 + line_coordinates_x[l] * M + k) = + -GeometryInfo<2>::d_linear_shape_function(px, l); + } + + { + constexpr std::array line_coordinates_z( + {{8, 9, 10, 11}}); + const Point<2> pz(p[0], p[1]); + for (unsigned int l = 0; l < 4; ++l) + lohvs(index_table, 8 + line_coordinates_z[l] * M + i) = + -GeometryInfo<2>::d_linear_shape_function(pz, l); + } + + // quads + lohvs(index_table, 8 + 12 * M + 0 * M * M + i * M + j) = + 1. - p[0]; + lohvs(index_table, 8 + 12 * M + 1 * M * M + i * M + j) = p[0]; + lohvs(index_table, 8 + 12 * M + 2 * M * M + k * M + i) = + 1. - p[1]; + lohvs(index_table, 8 + 12 * M + 3 * M * M + k * M + i) = p[1]; + lohvs(index_table, 8 + 12 * M + 4 * M * M + j * M + k) = + 1. - p[2]; + lohvs(index_table, 8 + 12 * M + 5 * M * M + j * M + k) = p[2]; + } + + // the sum of weights of the points at the outer rim should be one. + // check this + for (unsigned int unit_point = 0; unit_point < n_inner; ++unit_point) + Assert(std::fabs(std::accumulate(lohvs[unit_point].begin(), + lohvs[unit_point].end(), + 0.) - + 1) < 1e-13 * polynomial_degree, + ExcInternalError()); + + return lohvs; + } + + + + /** + * This function collects the output of + * compute_support_point_weights_on_{quad,hex} in a single data structure. + */ + inline std::vector> + compute_support_point_weights_perimeter_to_interior( + const unsigned int polynomial_degree, + const unsigned int dim) + { + Assert(dim > 0 && dim <= 3, ExcImpossibleInDim(dim)); + std::vector> output(dim); + if (polynomial_degree <= 1) + return output; + + // fill the 1D interior weights + QGaussLobatto<1> quadrature(polynomial_degree + 1); + output[0].reinit(polynomial_degree - 1, + GeometryInfo<1>::vertices_per_cell); + for (unsigned int q = 0; q < polynomial_degree - 1; ++q) + for (const unsigned int i : GeometryInfo<1>::vertex_indices()) + output[0](q, i) = + GeometryInfo<1>::d_linear_shape_function(quadrature.point(q + 1), + i); + + if (dim > 1) + output[1] = compute_support_point_weights_on_quad(polynomial_degree); + + if (dim > 2) + output[2] = compute_support_point_weights_on_hex(polynomial_degree); + + return output; + } + + + + /** + * Collects all interior points for the various dimensions. + */ + template + inline dealii::Table<2, double> + compute_support_point_weights_cell(const unsigned int polynomial_degree) + { + Assert(dim > 0 && dim <= 3, ExcImpossibleInDim(dim)); + if (polynomial_degree <= 1) + return dealii::Table<2, double>(); + + QGaussLobatto quadrature(polynomial_degree + 1); + const std::vector h2l = + FETools::hierarchic_to_lexicographic_numbering(polynomial_degree); + + dealii::Table<2, double> output(quadrature.size() - + GeometryInfo::vertices_per_cell, + GeometryInfo::vertices_per_cell); + for (unsigned int q = 0; q < output.size(0); ++q) + for (const unsigned int i : GeometryInfo::vertex_indices()) + output(q, i) = GeometryInfo::d_linear_shape_function( + quadrature.point(h2l[q + GeometryInfo::vertices_per_cell]), i); + + return output; + } + + + + /** + * Using the relative weights of the shape functions evaluated at + * one point on the reference cell (and stored in data.shape_values + * and accessed via data.shape(0,i)) and the locations of mapping + * support points (stored in data.mapping_support_points), compute + * the mapped location of that point in real space. + */ + template + inline Point + compute_mapped_location_of_point( + const typename dealii::MappingQGeneric::InternalData &data) + { + AssertDimension(data.shape_values.size(), + data.mapping_support_points.size()); + + // use now the InternalData to compute the point in real space. + Point p_real; + for (unsigned int i = 0; i < data.mapping_support_points.size(); ++i) + p_real += data.mapping_support_points[i] * data.shape(0, i); + + return p_real; + } + + + + /** + * Implementation of transform_real_to_unit_cell for either type double + * or VectorizedArray + */ + template + inline Point + do_transform_real_to_unit_cell_internal( + const Point & p, + const Point & initial_p_unit, + const std::vector> & points, + const std::vector> &polynomials_1d, + const std::vector & renumber, + const bool print_iterations_to_deallog = false) + { + AssertDimension(points.size(), + Utilities::pow(polynomials_1d.size(), dim)); + + // Newton iteration to solve + // f(x)=p(x)-p=0 + // where we are looking for 'x' and p(x) is the forward transformation + // from unit to real cell. We solve this using a Newton iteration + // x_{n+1}=x_n-[f'(x)]^{-1}f(x) + // The start value is set to be the linear approximation to the cell + + // The shape values and derivatives of the mapping at this point are + // previously computed. + + Point p_unit = initial_p_unit; + auto p_real = internal::evaluate_tensor_product_value_and_gradient( + polynomials_1d, points, p_unit, polynomials_1d.size() == 2, renumber); + + Tensor<1, spacedim, Number> f = p_real.first - p; + + // early out if we already have our point in all SIMD lanes, i.e., + // f.norm_square() < 1e-24 * p_real.second[0].norm_square(). To enable + // this step also for VectorizedArray where we do not have operator <, + // compare the result to zero which is defined for SIMD types + if (std::max(Number(0.), + f.norm_square() - 1e-24 * p_real.second[0].norm_square()) == + Number(0.)) + return p_unit; + + // we need to compare the position of the computed p(x) against the + // given point 'p'. We will terminate the iteration and return 'x' if + // they are less than eps apart. The question is how to choose eps -- + // or, put maybe more generally: in which norm we want these 'p' and + // 'p(x)' to be eps apart. + // + // the question is difficult since we may have to deal with very + // elongated cells where we may achieve 1e-12*h for the distance of + // these two points in the 'long' direction, but achieving this + // tolerance in the 'short' direction of the cell may not be possible + // + // what we do instead is then to terminate iterations if + // \| p(x) - p \|_A < eps + // where the A-norm is somehow induced by the transformation of the + // cell. in particular, we want to measure distances relative to the + // sizes of the cell in its principal directions. + // + // to define what exactly A should be, note that to first order we have + // the following (assuming that x* is the solution of the problem, i.e., + // p(x*)=p): + // p(x) - p = p(x) - p(x*) + // = -grad p(x) * (x*-x) + higher order terms + // This suggest to measure with a norm that corresponds to + // A = {[grad p(x]^T [grad p(x)]}^{-1} + // because then + // \| p(x) - p \|_A \approx \| x - x* \| + // Consequently, we will try to enforce that + // \| p(x) - p \|_A = \| f \| <= eps + // + // Note that using this norm is a bit dangerous since the norm changes + // in every iteration (A isn't fixed by depending on xk). However, if + // the cell is not too deformed (it may be stretched, but not twisted) + // then the mapping is almost linear and A is indeed constant or + // nearly so. + const double eps = 1.e-11; + const unsigned int newton_iteration_limit = 20; + + Point invalid_point; + invalid_point[0] = std::numeric_limits::infinity(); + bool try_project_to_unit_cell = false; + + unsigned int newton_iteration = 0; + Number f_weighted_norm_square = 1.; + Number last_f_weighted_norm_square = 1.; + + do + { + if (print_iterations_to_deallog) + deallog << "Newton iteration " << newton_iteration + << " for unit point guess " << p_unit << std::endl; + + // f'(x) + Tensor<2, spacedim, Number> df; + for (unsigned int d = 0; d < spacedim; ++d) + for (unsigned int e = 0; e < dim; ++e) + df[d][e] = p_real.second[e][d]; + + // check determinand(df) > 0 on all SIMD lanes + if (!(std::min(determinant(df), + Number(std::numeric_limits::min())) == + Number(std::numeric_limits::min()))) + { + // We allow to enter this function with an initial guess + // outside the unit cell. We might have run into invalid + // Jacobians due to that, so we should at least try once to go + // back to the unit cell and go on with the Newton iteration + // from there. Since the outside case is unlikely, we can + // afford spending the extra effort at this place. + if (try_project_to_unit_cell == false) + { + p_unit = GeometryInfo::project_to_unit_cell(p_unit); + p_real = internal::evaluate_tensor_product_value_and_gradient( + polynomials_1d, + points, + p_unit, + polynomials_1d.size() == 2, + renumber); + f = p_real.first - p; + f_weighted_norm_square = 1.; + last_f_weighted_norm_square = 1; + try_project_to_unit_cell = true; + continue; + } + else + return invalid_point; + } + + // Solve [f'(x)]d=f(x) + const Tensor<2, spacedim, Number> df_inverse = invert(df); + const Tensor<1, spacedim, Number> delta = df_inverse * f; + last_f_weighted_norm_square = delta.norm_square(); + + if (print_iterations_to_deallog) + deallog << " delta=" << delta << std::endl; + + // do a line search + double step_length = 1; + do + { + // update of p_unit. The spacedim-th component of transformed + // point is simply ignored in codimension one case. When this + // component is not zero, then we are projecting the point to + // the surface or curve identified by the cell. + Point p_unit_trial = p_unit; + for (unsigned int i = 0; i < dim; ++i) + p_unit_trial[i] -= step_length * delta[i]; + + // shape values and derivatives at new p_unit point + const auto p_real_trial = + internal::evaluate_tensor_product_value_and_gradient( + polynomials_1d, + points, + p_unit_trial, + polynomials_1d.size() == 2, + renumber); + const Tensor<1, spacedim, Number> f_trial = + p_real_trial.first - p; + f_weighted_norm_square = (df_inverse * f_trial).norm_square(); + + if (print_iterations_to_deallog) + deallog << " step_length=" << step_length << std::endl + << " ||f || =" << f.norm() << std::endl + << " ||f*|| =" << f_trial.norm() << std::endl + << " ||f*||_A =" + << std::sqrt(f_weighted_norm_square) << std::endl; + + // See if we are making progress with the current step length + // and if not, reduce it by a factor of two and try again. + // + // Strictly speaking, we should probably use the same norm as we + // use for the outer algorithm. In practice, line search is just + // a crutch to find a "reasonable" step length, and so using the + // l2 norm is probably just fine. + // + // check f_trial.norm() < f.norm() in SIMD form. This is a bit + // more complicated because some SIMD lanes might not be doing + // any progress any more as they have already reached roundoff + // accuracy. We define that as the case of updates less than + // 1e-6. The tolerance might seem coarse but since we are + // dealing with a Newton iteration of a polynomial function we + // either converge quadratically or not any more. Thus, our + // selection is to terminate if either the norm of f is + // decreasing or that threshold of 1e-6 is reached. + if (std::max(f_weighted_norm_square - 1e-6 * 1e-6, Number(0.)) * + std::max(f_trial.norm_square() - f.norm_square(), + Number(0.)) == + Number(0.)) + { + p_real = p_real_trial; + p_unit = p_unit_trial; + f = f_trial; + break; + } + else if (step_length > 0.05) + step_length *= 0.5; + else + break; + } + while (true); + + // In case we terminated the line search due to the step becoming + // too small, we give the iteration another try with the + // projection of the initial guess to the unit cell before we give + // up, just like for the negative determinant case. + if (step_length <= 0.05 && try_project_to_unit_cell == false) + { + p_unit = GeometryInfo::project_to_unit_cell(p_unit); + p_real = internal::evaluate_tensor_product_value_and_gradient( + polynomials_1d, + points, + p_unit, + polynomials_1d.size() == 2, + renumber); + f = p_real.first - p; + f_weighted_norm_square = 1.; + last_f_weighted_norm_square = 1; + try_project_to_unit_cell = true; + continue; + } + else if (step_length <= 0.05) + return invalid_point; + + ++newton_iteration; + if (newton_iteration > newton_iteration_limit) + return invalid_point; + } + // Stop if f_weighted_norm_square <= eps^2 on all SIMD lanes or if the + // weighted norm is less than 1e-6 and the improvement against the + // previous step was less than a factor of 10 (in that regime, we + // either have quadratic convergence or roundoff errors due to a bad + // mapping) + while ( + !(std::max(f_weighted_norm_square - eps * eps, Number(0.)) * + std::max(last_f_weighted_norm_square - + std::min(f_weighted_norm_square, Number(1e-6 * 1e-6)) * + 100., + Number(0.)) == + Number(0.))); + + if (print_iterations_to_deallog) + deallog << "Iteration converged for p_unit = [ " << p_unit + << " ] and iteration error " + << std::sqrt(f_weighted_norm_square) << std::endl; + + return p_unit; + } + + + + /** + * Implementation of transform_real_to_unit_cell for dim==spacedim-1 + */ + template + inline Point + do_transform_real_to_unit_cell_internal_codim1( + const typename dealii::Triangulation::cell_iterator &cell, + const Point & p, + const Point &initial_p_unit, + typename dealii::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 < 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 &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(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 &>(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>(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 &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(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::ExcTransformationFailed())); + + return p_unit; + } + + /** + * In case the quadrature formula is a tensor product, this is a + * replacement for maybe_compute_q_points(), maybe_update_Jacobians() and + * maybe_update_jacobian_grads() + */ + template + inline void + maybe_update_q_points_Jacobians_and_grads_tensor( + const CellSimilarity::Similarity cell_similarity, + const typename dealii::MappingQGeneric::InternalData &data, + std::vector> & quadrature_points, + std::vector> &jacobian_grads) + { + const UpdateFlags update_flags = data.update_each; + + const unsigned int n_shape_values = data.n_shape_functions; + const unsigned int n_q_points = data.shape_info.n_q_points; + constexpr unsigned int n_lanes = VectorizedArray::size(); + constexpr unsigned int n_comp = 1 + (spacedim - 1) / n_lanes; + constexpr unsigned int n_hessians = (dim * (dim + 1)) / 2; + + EvaluationFlags::EvaluationFlags evaluation_flag = + (update_flags & update_quadrature_points ? EvaluationFlags::values : + EvaluationFlags::nothing) | + ((cell_similarity != CellSimilarity::translation) && + (update_flags & update_contravariant_transformation) ? + EvaluationFlags::gradients : + EvaluationFlags::nothing) | + ((cell_similarity != CellSimilarity::translation) && + (update_flags & update_jacobian_grads) ? + EvaluationFlags::hessians : + EvaluationFlags::nothing); + + Assert(!(evaluation_flag & EvaluationFlags::values) || n_q_points > 0, + ExcInternalError()); + Assert(!(evaluation_flag & EvaluationFlags::values) || + n_q_points == quadrature_points.size(), + ExcDimensionMismatch(n_q_points, quadrature_points.size())); + Assert(!(evaluation_flag & EvaluationFlags::gradients) || + data.n_shape_functions > 0, + ExcInternalError()); + Assert(!(evaluation_flag & EvaluationFlags::gradients) || + n_q_points == data.contravariant.size(), + ExcDimensionMismatch(n_q_points, data.contravariant.size())); + Assert(!(evaluation_flag & EvaluationFlags::hessians) || + n_q_points == jacobian_grads.size(), + ExcDimensionMismatch(n_q_points, jacobian_grads.size())); + + // shortcut in case we have an identity interpolation and only request + // the quadrature points + if (evaluation_flag == EvaluationFlags::values && + data.shape_info.element_type == + internal::MatrixFreeFunctions::tensor_symmetric_collocation) + { + for (unsigned int q = 0; q < n_q_points; ++q) + quadrature_points[q] = + data.mapping_support_points[data.shape_info + .lexicographic_numbering[q]]; + return; + } + + // prepare arrays + if (evaluation_flag != EvaluationFlags::nothing) + { + data.values_dofs.resize(n_comp * n_shape_values); + data.values_quad.resize(n_comp * n_q_points); + data.gradients_quad.resize(n_comp * n_q_points * dim); + data.scratch.resize(2 * std::max(n_q_points, n_shape_values)); + + if (evaluation_flag & EvaluationFlags::hessians) + data.hessians_quad.resize(n_comp * n_q_points * n_hessians); + + const std::vector &renumber_to_lexicographic = + data.shape_info.lexicographic_numbering; + for (unsigned int i = 0; i < n_shape_values; ++i) + for (unsigned int d = 0; d < spacedim; ++d) + { + const unsigned int in_comp = d % n_lanes; + const unsigned int out_comp = d / n_lanes; + data.values_dofs[out_comp * n_shape_values + i][in_comp] = + data.mapping_support_points[renumber_to_lexicographic[i]][d]; + } + + // do the actual tensorized evaluation + internal::FEEvaluationFactory>:: + evaluate(n_comp, + evaluation_flag, + data.shape_info, + data.values_dofs.begin(), + data.values_quad.begin(), + data.gradients_quad.begin(), + data.hessians_quad.begin(), + data.scratch.begin()); + } + + // do the postprocessing + if (evaluation_flag & EvaluationFlags::values) + { + for (unsigned int out_comp = 0; out_comp < n_comp; ++out_comp) + for (unsigned int i = 0; i < n_q_points; ++i) + for (unsigned int in_comp = 0; + in_comp < n_lanes && in_comp < spacedim - out_comp * n_lanes; + ++in_comp) + quadrature_points[i][out_comp * n_lanes + in_comp] = + data.values_quad[out_comp * n_q_points + i][in_comp]; + } + + if (evaluation_flag & EvaluationFlags::gradients) + { + std::fill(data.contravariant.begin(), + data.contravariant.end(), + DerivativeForm<1, dim, spacedim>()); + // We need to reinterpret the data after evaluate has been applied. + for (unsigned int out_comp = 0; out_comp < n_comp; ++out_comp) + for (unsigned int point = 0; point < n_q_points; ++point) + for (unsigned int j = 0; j < dim; ++j) + for (unsigned int in_comp = 0; + in_comp < n_lanes && + in_comp < spacedim - out_comp * n_lanes; + ++in_comp) + { + const unsigned int total_number = point * dim + j; + const unsigned int new_comp = total_number / n_q_points; + const unsigned int new_point = total_number % n_q_points; + data.contravariant[new_point][out_comp * n_lanes + in_comp] + [new_comp] = + data + .gradients_quad[(out_comp * n_q_points + point) * dim + + j][in_comp]; + } + } + if (update_flags & update_covariant_transformation) + if (cell_similarity != CellSimilarity::translation) + for (unsigned int point = 0; point < n_q_points; ++point) + data.covariant[point] = + (data.contravariant[point]).covariant_form(); + + if (update_flags & update_volume_elements) + if (cell_similarity != CellSimilarity::translation) + for (unsigned int point = 0; point < n_q_points; ++point) + data.volume_elements[point] = + data.contravariant[point].determinant(); + + if (evaluation_flag & EvaluationFlags::hessians) + { + constexpr int desymmetrize_3d[6][2] = { + {0, 0}, {1, 1}, {2, 2}, {0, 1}, {0, 2}, {1, 2}}; + constexpr int desymmetrize_2d[3][2] = {{0, 0}, {1, 1}, {0, 1}}; + + // We need to reinterpret the data after evaluate has been applied. + for (unsigned int out_comp = 0; out_comp < n_comp; ++out_comp) + for (unsigned int point = 0; point < n_q_points; ++point) + for (unsigned int j = 0; j < n_hessians; ++j) + for (unsigned int in_comp = 0; + in_comp < n_lanes && + in_comp < spacedim - out_comp * n_lanes; + ++in_comp) + { + const unsigned int total_number = point * n_hessians + j; + const unsigned int new_point = total_number % n_q_points; + const unsigned int new_hessian_comp = + total_number / n_q_points; + const unsigned int new_hessian_comp_i = + dim == 2 ? desymmetrize_2d[new_hessian_comp][0] : + desymmetrize_3d[new_hessian_comp][0]; + const unsigned int new_hessian_comp_j = + dim == 2 ? desymmetrize_2d[new_hessian_comp][1] : + desymmetrize_3d[new_hessian_comp][1]; + const double value = + data.hessians_quad[(out_comp * n_q_points + point) * + n_hessians + + j][in_comp]; + jacobian_grads[new_point][out_comp * n_lanes + in_comp] + [new_hessian_comp_i][new_hessian_comp_j] = + value; + jacobian_grads[new_point][out_comp * n_lanes + in_comp] + [new_hessian_comp_j][new_hessian_comp_i] = + value; + } + } + } + + + /** + * Compute the locations of quadrature points on the object described by + * the first argument (and the cell for which the mapping support points + * have already been set), but only if the update_flags of the @p data + * argument indicate so. + */ + template + inline void + maybe_compute_q_points( + const typename QProjector::DataSetDescriptor data_set, + const typename dealii::MappingQGeneric::InternalData &data, + std::vector> &quadrature_points) + { + const UpdateFlags update_flags = data.update_each; + + if (update_flags & update_quadrature_points) + for (unsigned int point = 0; point < quadrature_points.size(); ++point) + { + const double * shape = &data.shape(point + data_set, 0); + Point result = + (shape[0] * data.mapping_support_points[0]); + for (unsigned int k = 1; k < data.n_shape_functions; ++k) + for (unsigned int i = 0; i < spacedim; ++i) + result[i] += shape[k] * data.mapping_support_points[k][i]; + quadrature_points[point] = result; + } + } + + + + /** + * Update the co- and contravariant matrices as well as their determinant, + * for the cell + * described stored in the data object, but only if the update_flags of the @p data + * argument indicate so. + * + * Skip the computation if possible as indicated by the first argument. + */ + template + inline void + maybe_update_Jacobians( + const CellSimilarity::Similarity cell_similarity, + const typename dealii::QProjector::DataSetDescriptor data_set, + const typename dealii::MappingQGeneric::InternalData &data) + { + const UpdateFlags update_flags = data.update_each; + + if (update_flags & update_contravariant_transformation) + // if the current cell is just a + // translation of the previous one, no + // need to recompute jacobians... + if (cell_similarity != CellSimilarity::translation) + { + const unsigned int n_q_points = data.contravariant.size(); + + std::fill(data.contravariant.begin(), + data.contravariant.end(), + DerivativeForm<1, dim, spacedim>()); + + Assert(data.n_shape_functions > 0, ExcInternalError()); + + const Tensor<1, spacedim> *supp_pts = + data.mapping_support_points.data(); + + for (unsigned int point = 0; point < n_q_points; ++point) + { + const Tensor<1, dim> *data_derv = + &data.derivative(point + data_set, 0); + + double result[spacedim][dim]; + + // peel away part of sum to avoid zeroing the + // entries and adding for the first time + for (unsigned int i = 0; i < spacedim; ++i) + for (unsigned int j = 0; j < dim; ++j) + result[i][j] = data_derv[0][j] * supp_pts[0][i]; + for (unsigned int k = 1; k < data.n_shape_functions; ++k) + for (unsigned int i = 0; i < spacedim; ++i) + for (unsigned int j = 0; j < dim; ++j) + result[i][j] += data_derv[k][j] * supp_pts[k][i]; + + // write result into contravariant data. for + // j=dim in the case dim + inline void + maybe_update_jacobian_grads( + const CellSimilarity::Similarity cell_similarity, + const typename QProjector::DataSetDescriptor data_set, + const typename dealii::MappingQGeneric::InternalData &data, + std::vector> &jacobian_grads) + { + const UpdateFlags update_flags = data.update_each; + if (update_flags & update_jacobian_grads) + { + const unsigned int n_q_points = jacobian_grads.size(); + + if (cell_similarity != CellSimilarity::translation) + for (unsigned int point = 0; point < n_q_points; ++point) + { + const Tensor<2, dim> *second = + &data.second_derivative(point + data_set, 0); + double result[spacedim][dim][dim]; + for (unsigned int i = 0; i < spacedim; ++i) + for (unsigned int j = 0; j < dim; ++j) + for (unsigned int l = 0; l < dim; ++l) + result[i][j][l] = + (second[0][j][l] * data.mapping_support_points[0][i]); + for (unsigned int k = 1; k < data.n_shape_functions; ++k) + for (unsigned int i = 0; i < spacedim; ++i) + for (unsigned int j = 0; j < dim; ++j) + for (unsigned int l = 0; l < dim; ++l) + result[i][j][l] += + (second[k][j][l] * data.mapping_support_points[k][i]); + + for (unsigned int i = 0; i < spacedim; ++i) + for (unsigned int j = 0; j < dim; ++j) + for (unsigned int l = 0; l < dim; ++l) + jacobian_grads[point][i][j][l] = result[i][j][l]; + } + } + } + + + + /** + * Update the Hessian of the transformation from unit to real cell, the + * Jacobian gradients, pushed forward to the real cell coordinates. + * + * Skip the computation if possible as indicated by the first argument. + */ + template + inline void + maybe_update_jacobian_pushed_forward_grads( + const CellSimilarity::Similarity cell_similarity, + const typename QProjector::DataSetDescriptor data_set, + const typename dealii::MappingQGeneric::InternalData &data, + std::vector> &jacobian_pushed_forward_grads) + { + const UpdateFlags update_flags = data.update_each; + if (update_flags & update_jacobian_pushed_forward_grads) + { + const unsigned int n_q_points = jacobian_pushed_forward_grads.size(); + + if (cell_similarity != CellSimilarity::translation) + { + double tmp[spacedim][spacedim][spacedim]; + for (unsigned int point = 0; point < n_q_points; ++point) + { + const Tensor<2, dim> *second = + &data.second_derivative(point + data_set, 0); + double result[spacedim][dim][dim]; + for (unsigned int i = 0; i < spacedim; ++i) + for (unsigned int j = 0; j < dim; ++j) + for (unsigned int l = 0; l < dim; ++l) + result[i][j][l] = + (second[0][j][l] * data.mapping_support_points[0][i]); + for (unsigned int k = 1; k < data.n_shape_functions; ++k) + for (unsigned int i = 0; i < spacedim; ++i) + for (unsigned int j = 0; j < dim; ++j) + for (unsigned int l = 0; l < dim; ++l) + result[i][j][l] += + (second[k][j][l] * + data.mapping_support_points[k][i]); + + // first push forward the j-components + for (unsigned int i = 0; i < spacedim; ++i) + for (unsigned int j = 0; j < spacedim; ++j) + for (unsigned int l = 0; l < dim; ++l) + { + tmp[i][j][l] = + result[i][0][l] * data.covariant[point][j][0]; + for (unsigned int jr = 1; jr < dim; ++jr) + { + tmp[i][j][l] += + result[i][jr][l] * data.covariant[point][j][jr]; + } + } + + // now, pushing forward the l-components + for (unsigned int i = 0; i < spacedim; ++i) + for (unsigned int j = 0; j < spacedim; ++j) + for (unsigned int l = 0; l < spacedim; ++l) + { + jacobian_pushed_forward_grads[point][i][j][l] = + tmp[i][j][0] * data.covariant[point][l][0]; + for (unsigned int lr = 1; lr < dim; ++lr) + { + jacobian_pushed_forward_grads[point][i][j][l] += + tmp[i][j][lr] * data.covariant[point][l][lr]; + } + } + } + } + } + } + + + + /** + * Update the third derivatives of the transformation from unit to real + * cell, the Jacobian hessians. + * + * Skip the computation if possible as indicated by the first argument. + */ + template + inline void + maybe_update_jacobian_2nd_derivatives( + const CellSimilarity::Similarity cell_similarity, + const typename QProjector::DataSetDescriptor data_set, + const typename dealii::MappingQGeneric::InternalData &data, + std::vector> &jacobian_2nd_derivatives) + { + const UpdateFlags update_flags = data.update_each; + if (update_flags & update_jacobian_2nd_derivatives) + { + const unsigned int n_q_points = jacobian_2nd_derivatives.size(); + + if (cell_similarity != CellSimilarity::translation) + { + for (unsigned int point = 0; point < n_q_points; ++point) + { + const Tensor<3, dim> *third = + &data.third_derivative(point + data_set, 0); + double result[spacedim][dim][dim][dim]; + for (unsigned int i = 0; i < spacedim; ++i) + for (unsigned int j = 0; j < dim; ++j) + for (unsigned int l = 0; l < dim; ++l) + for (unsigned int m = 0; m < dim; ++m) + result[i][j][l][m] = + (third[0][j][l][m] * + data.mapping_support_points[0][i]); + for (unsigned int k = 1; k < data.n_shape_functions; ++k) + for (unsigned int i = 0; i < spacedim; ++i) + for (unsigned int j = 0; j < dim; ++j) + for (unsigned int l = 0; l < dim; ++l) + for (unsigned int m = 0; m < dim; ++m) + result[i][j][l][m] += + (third[k][j][l][m] * + data.mapping_support_points[k][i]); + + for (unsigned int i = 0; i < spacedim; ++i) + for (unsigned int j = 0; j < dim; ++j) + for (unsigned int l = 0; l < dim; ++l) + for (unsigned int m = 0; m < dim; ++m) + jacobian_2nd_derivatives[point][i][j][l][m] = + result[i][j][l][m]; + } + } + } + } + + + + /** + * Update the Hessian of the Hessian of the transformation from unit + * to real cell, the Jacobian Hessian gradients, pushed forward to the + * real cell coordinates. + * + * Skip the computation if possible as indicated by the first argument. + */ + template + inline void + maybe_update_jacobian_pushed_forward_2nd_derivatives( + const CellSimilarity::Similarity cell_similarity, + const typename QProjector::DataSetDescriptor data_set, + const typename dealii::MappingQGeneric::InternalData &data, + std::vector> &jacobian_pushed_forward_2nd_derivatives) + { + const UpdateFlags update_flags = data.update_each; + if (update_flags & update_jacobian_pushed_forward_2nd_derivatives) + { + const unsigned int n_q_points = + jacobian_pushed_forward_2nd_derivatives.size(); + + if (cell_similarity != CellSimilarity::translation) + { + double tmp[spacedim][spacedim][spacedim][spacedim]; + for (unsigned int point = 0; point < n_q_points; ++point) + { + const Tensor<3, dim> *third = + &data.third_derivative(point + data_set, 0); + double result[spacedim][dim][dim][dim]; + for (unsigned int i = 0; i < spacedim; ++i) + for (unsigned int j = 0; j < dim; ++j) + for (unsigned int l = 0; l < dim; ++l) + for (unsigned int m = 0; m < dim; ++m) + result[i][j][l][m] = + (third[0][j][l][m] * + data.mapping_support_points[0][i]); + for (unsigned int k = 1; k < data.n_shape_functions; ++k) + for (unsigned int i = 0; i < spacedim; ++i) + for (unsigned int j = 0; j < dim; ++j) + for (unsigned int l = 0; l < dim; ++l) + for (unsigned int m = 0; m < dim; ++m) + result[i][j][l][m] += + (third[k][j][l][m] * + data.mapping_support_points[k][i]); + + // push forward the j-coordinate + for (unsigned int i = 0; i < spacedim; ++i) + for (unsigned int j = 0; j < spacedim; ++j) + for (unsigned int l = 0; l < dim; ++l) + for (unsigned int m = 0; m < dim; ++m) + { + jacobian_pushed_forward_2nd_derivatives + [point][i][j][l][m] = result[i][0][l][m] * + data.covariant[point][j][0]; + for (unsigned int jr = 1; jr < dim; ++jr) + jacobian_pushed_forward_2nd_derivatives[point][i] + [j][l] + [m] += + result[i][jr][l][m] * + data.covariant[point][j][jr]; + } + + // push forward the l-coordinate + for (unsigned int i = 0; i < spacedim; ++i) + for (unsigned int j = 0; j < spacedim; ++j) + for (unsigned int l = 0; l < spacedim; ++l) + for (unsigned int m = 0; m < dim; ++m) + { + tmp[i][j][l][m] = + jacobian_pushed_forward_2nd_derivatives[point][i] + [j][0][m] * + data.covariant[point][l][0]; + for (unsigned int lr = 1; lr < dim; ++lr) + tmp[i][j][l][m] += + jacobian_pushed_forward_2nd_derivatives[point] + [i][j] + [lr][m] * + data.covariant[point][l][lr]; + } + + // push forward the m-coordinate + for (unsigned int i = 0; i < spacedim; ++i) + for (unsigned int j = 0; j < spacedim; ++j) + for (unsigned int l = 0; l < spacedim; ++l) + for (unsigned int m = 0; m < spacedim; ++m) + { + jacobian_pushed_forward_2nd_derivatives + [point][i][j][l][m] = + tmp[i][j][l][0] * data.covariant[point][m][0]; + for (unsigned int mr = 1; mr < dim; ++mr) + jacobian_pushed_forward_2nd_derivatives[point][i] + [j][l] + [m] += + tmp[i][j][l][mr] * data.covariant[point][m][mr]; + } + } + } + } + } + + + + /** + * Update the fourth derivatives of the transformation from unit to real + * cell, the Jacobian hessian gradients. + * + * Skip the computation if possible as indicated by the first argument. + */ + template + inline void + maybe_update_jacobian_3rd_derivatives( + const CellSimilarity::Similarity cell_similarity, + const typename QProjector::DataSetDescriptor data_set, + const typename dealii::MappingQGeneric::InternalData &data, + std::vector> &jacobian_3rd_derivatives) + { + const UpdateFlags update_flags = data.update_each; + if (update_flags & update_jacobian_3rd_derivatives) + { + const unsigned int n_q_points = jacobian_3rd_derivatives.size(); + + if (cell_similarity != CellSimilarity::translation) + { + for (unsigned int point = 0; point < n_q_points; ++point) + { + const Tensor<4, dim> *fourth = + &data.fourth_derivative(point + data_set, 0); + double result[spacedim][dim][dim][dim][dim]; + for (unsigned int i = 0; i < spacedim; ++i) + for (unsigned int j = 0; j < dim; ++j) + for (unsigned int l = 0; l < dim; ++l) + for (unsigned int m = 0; m < dim; ++m) + for (unsigned int n = 0; n < dim; ++n) + result[i][j][l][m][n] = + (fourth[0][j][l][m][n] * + data.mapping_support_points[0][i]); + for (unsigned int k = 1; k < data.n_shape_functions; ++k) + for (unsigned int i = 0; i < spacedim; ++i) + for (unsigned int j = 0; j < dim; ++j) + for (unsigned int l = 0; l < dim; ++l) + for (unsigned int m = 0; m < dim; ++m) + for (unsigned int n = 0; n < dim; ++n) + result[i][j][l][m][n] += + (fourth[k][j][l][m][n] * + data.mapping_support_points[k][i]); + + for (unsigned int i = 0; i < spacedim; ++i) + for (unsigned int j = 0; j < dim; ++j) + for (unsigned int l = 0; l < dim; ++l) + for (unsigned int m = 0; m < dim; ++m) + for (unsigned int n = 0; n < dim; ++n) + jacobian_3rd_derivatives[point][i][j][l][m][n] = + result[i][j][l][m][n]; + } + } + } + } + + + + /** + * Update the Hessian gradient of the transformation from unit to real + * cell, the Jacobian Hessians, pushed forward to the real cell + * coordinates. + * + * Skip the computation if possible as indicated by the first argument. + */ + template + inline void + maybe_update_jacobian_pushed_forward_3rd_derivatives( + const CellSimilarity::Similarity cell_similarity, + const typename QProjector::DataSetDescriptor data_set, + const typename dealii::MappingQGeneric::InternalData &data, + std::vector> &jacobian_pushed_forward_3rd_derivatives) + { + const UpdateFlags update_flags = data.update_each; + if (update_flags & update_jacobian_pushed_forward_3rd_derivatives) + { + const unsigned int n_q_points = + jacobian_pushed_forward_3rd_derivatives.size(); + + if (cell_similarity != CellSimilarity::translation) + { + double tmp[spacedim][spacedim][spacedim][spacedim][spacedim]; + for (unsigned int point = 0; point < n_q_points; ++point) + { + const Tensor<4, dim> *fourth = + &data.fourth_derivative(point + data_set, 0); + double result[spacedim][dim][dim][dim][dim]; + for (unsigned int i = 0; i < spacedim; ++i) + for (unsigned int j = 0; j < dim; ++j) + for (unsigned int l = 0; l < dim; ++l) + for (unsigned int m = 0; m < dim; ++m) + for (unsigned int n = 0; n < dim; ++n) + result[i][j][l][m][n] = + (fourth[0][j][l][m][n] * + data.mapping_support_points[0][i]); + for (unsigned int k = 1; k < data.n_shape_functions; ++k) + for (unsigned int i = 0; i < spacedim; ++i) + for (unsigned int j = 0; j < dim; ++j) + for (unsigned int l = 0; l < dim; ++l) + for (unsigned int m = 0; m < dim; ++m) + for (unsigned int n = 0; n < dim; ++n) + result[i][j][l][m][n] += + (fourth[k][j][l][m][n] * + data.mapping_support_points[k][i]); + + // push-forward the j-coordinate + for (unsigned int i = 0; i < spacedim; ++i) + for (unsigned int j = 0; j < spacedim; ++j) + for (unsigned int l = 0; l < dim; ++l) + for (unsigned int m = 0; m < dim; ++m) + for (unsigned int n = 0; n < dim; ++n) + { + tmp[i][j][l][m][n] = result[i][0][l][m][n] * + data.covariant[point][j][0]; + for (unsigned int jr = 1; jr < dim; ++jr) + tmp[i][j][l][m][n] += + result[i][jr][l][m][n] * + data.covariant[point][j][jr]; + } + + // push-forward the l-coordinate + for (unsigned int i = 0; i < spacedim; ++i) + for (unsigned int j = 0; j < spacedim; ++j) + for (unsigned int l = 0; l < spacedim; ++l) + for (unsigned int m = 0; m < dim; ++m) + for (unsigned int n = 0; n < dim; ++n) + { + jacobian_pushed_forward_3rd_derivatives + [point][i][j][l][m][n] = + tmp[i][j][0][m][n] * + data.covariant[point][l][0]; + for (unsigned int lr = 1; lr < dim; ++lr) + jacobian_pushed_forward_3rd_derivatives[point] + [i][j][l] + [m][n] += + tmp[i][j][lr][m][n] * + data.covariant[point][l][lr]; + } + + // push-forward the m-coordinate + for (unsigned int i = 0; i < spacedim; ++i) + for (unsigned int j = 0; j < spacedim; ++j) + for (unsigned int l = 0; l < spacedim; ++l) + for (unsigned int m = 0; m < spacedim; ++m) + for (unsigned int n = 0; n < dim; ++n) + { + tmp[i][j][l][m][n] = + jacobian_pushed_forward_3rd_derivatives[point] + [i][j][l] + [0][n] * + data.covariant[point][m][0]; + for (unsigned int mr = 1; mr < dim; ++mr) + tmp[i][j][l][m][n] += + jacobian_pushed_forward_3rd_derivatives + [point][i][j][l][mr][n] * + data.covariant[point][m][mr]; + } + + // push-forward the n-coordinate + for (unsigned int i = 0; i < spacedim; ++i) + for (unsigned int j = 0; j < spacedim; ++j) + for (unsigned int l = 0; l < spacedim; ++l) + for (unsigned int m = 0; m < spacedim; ++m) + for (unsigned int n = 0; n < spacedim; ++n) + { + jacobian_pushed_forward_3rd_derivatives + [point][i][j][l][m][n] = + tmp[i][j][l][m][0] * + data.covariant[point][n][0]; + for (unsigned int nr = 1; nr < dim; ++nr) + jacobian_pushed_forward_3rd_derivatives[point] + [i][j][l] + [m][n] += + tmp[i][j][l][m][nr] * + data.covariant[point][n][nr]; + } + } + } + } + } + + + + /** + * Depending on what information is called for in the update flags of the + * @p data object, compute the various pieces of information that is + * required by the fill_fe_face_values() and fill_fe_subface_values() + * functions. This function simply unifies the work that would be done by + * those two functions. + * + * The resulting data is put into the @p output_data argument. + */ + template + inline void + maybe_compute_face_data( + const dealii::MappingQGeneric &mapping, + const typename dealii::Triangulation::cell_iterator &cell, + const unsigned int face_no, + const unsigned int subface_no, + const unsigned int n_q_points, + const std::vector &weights, + const typename dealii::MappingQGeneric::InternalData &data, + internal::FEValuesImplementation::MappingRelatedData + &output_data) + { + const UpdateFlags update_flags = data.update_each; + + if (update_flags & + (update_boundary_forms | update_normal_vectors | update_jacobians | + update_JxW_values | update_inverse_jacobians)) + { + if (update_flags & update_boundary_forms) + AssertDimension(output_data.boundary_forms.size(), n_q_points); + if (update_flags & update_normal_vectors) + AssertDimension(output_data.normal_vectors.size(), n_q_points); + if (update_flags & update_JxW_values) + AssertDimension(output_data.JxW_values.size(), n_q_points); + + Assert(data.aux.size() + 1 >= dim, ExcInternalError()); + + // first compute some common data that is used for evaluating + // all of the flags below + + // map the unit tangentials to the real cell. checking for d!=dim-1 + // eliminates compiler warnings regarding unsigned int expressions < + // 0. + for (unsigned int d = 0; d != dim - 1; ++d) + { + Assert(face_no + GeometryInfo::faces_per_cell * d < + data.unit_tangentials.size(), + ExcInternalError()); + Assert( + data.aux[d].size() <= + data + .unit_tangentials[face_no + + GeometryInfo::faces_per_cell * d] + .size(), + ExcInternalError()); + + mapping.transform( + make_array_view( + data.unit_tangentials[face_no + + GeometryInfo::faces_per_cell * d]), + mapping_contravariant, + data, + make_array_view(data.aux[d])); + } + + if (update_flags & update_boundary_forms) + { + // if dim==spacedim, we can use the unit tangentials to compute + // the boundary form by simply taking the cross product + if (dim == spacedim) + { + for (unsigned int i = 0; i < n_q_points; ++i) + switch (dim) + { + case 1: + // in 1d, we don't have access to any of the + // data.aux fields (because it has only dim-1 + // components), but we can still compute the + // boundary form by simply looking at the number of + // the face + output_data.boundary_forms[i][0] = + (face_no == 0 ? -1 : +1); + break; + case 2: + output_data.boundary_forms[i] = + cross_product_2d(data.aux[0][i]); + break; + case 3: + output_data.boundary_forms[i] = + cross_product_3d(data.aux[0][i], data.aux[1][i]); + break; + default: + Assert(false, ExcNotImplemented()); + } + } + else //(dim < spacedim) + { + // in the codim-one case, the boundary form results from the + // cross product of all the face tangential vectors and the + // cell normal vector + // + // to compute the cell normal, use the same method used in + // fill_fe_values for cells above + AssertDimension(data.contravariant.size(), n_q_points); + + for (unsigned int point = 0; point < n_q_points; ++point) + { + if (dim == 1) + { + // J is a tangent vector + output_data.boundary_forms[point] = + data.contravariant[point].transpose()[0]; + output_data.boundary_forms[point] /= + (face_no == 0 ? -1. : +1.) * + output_data.boundary_forms[point].norm(); + } + + if (dim == 2) + { + const DerivativeForm<1, spacedim, dim> DX_t = + data.contravariant[point].transpose(); + + Tensor<1, spacedim> cell_normal = + cross_product_3d(DX_t[0], DX_t[1]); + cell_normal /= cell_normal.norm(); + + // then compute the face normal from the face + // tangent and the cell normal: + output_data.boundary_forms[point] = + cross_product_3d(data.aux[0][point], cell_normal); + } + } + } + } + + if (update_flags & update_JxW_values) + for (unsigned int i = 0; i < output_data.boundary_forms.size(); ++i) + { + output_data.JxW_values[i] = + output_data.boundary_forms[i].norm() * weights[i]; + + if (subface_no != numbers::invalid_unsigned_int) + { + const double area_ratio = GeometryInfo::subface_ratio( + cell->subface_case(face_no), subface_no); + output_data.JxW_values[i] *= area_ratio; + } + } + + if (update_flags & update_normal_vectors) + for (unsigned int i = 0; i < output_data.normal_vectors.size(); ++i) + output_data.normal_vectors[i] = + Point(output_data.boundary_forms[i] / + output_data.boundary_forms[i].norm()); + + if (update_flags & update_jacobians) + for (unsigned int point = 0; point < n_q_points; ++point) + output_data.jacobians[point] = data.contravariant[point]; + + if (update_flags & update_inverse_jacobians) + for (unsigned int point = 0; point < n_q_points; ++point) + output_data.inverse_jacobians[point] = + data.covariant[point].transpose(); + } + } + + + /** + * Do the work of MappingQGeneric::fill_fe_face_values() and + * MappingQGeneric::fill_fe_subface_values() in a generic way, + * using the 'data_set' to differentiate whether we will + * work on a face (and if so, which one) or subface. + */ + template + inline void + do_fill_fe_face_values( + const dealii::MappingQGeneric &mapping, + const typename dealii::Triangulation::cell_iterator &cell, + const unsigned int face_no, + const unsigned int subface_no, + const typename QProjector::DataSetDescriptor data_set, + const Quadrature & quadrature, + const typename dealii::MappingQGeneric::InternalData &data, + internal::FEValuesImplementation::MappingRelatedData + &output_data) + { + if (dim > 1 && data.tensor_product_quadrature) + { + maybe_update_q_points_Jacobians_and_grads_tensor( + CellSimilarity::none, + data, + output_data.quadrature_points, + output_data.jacobian_grads); + } + else + { + maybe_compute_q_points(data_set, + data, + output_data.quadrature_points); + maybe_update_Jacobians(CellSimilarity::none, + data_set, + data); + maybe_update_jacobian_grads( + CellSimilarity::none, data_set, data, output_data.jacobian_grads); + } + maybe_update_jacobian_pushed_forward_grads( + CellSimilarity::none, + data_set, + data, + output_data.jacobian_pushed_forward_grads); + maybe_update_jacobian_2nd_derivatives( + CellSimilarity::none, + data_set, + data, + output_data.jacobian_2nd_derivatives); + maybe_update_jacobian_pushed_forward_2nd_derivatives( + CellSimilarity::none, + data_set, + data, + output_data.jacobian_pushed_forward_2nd_derivatives); + maybe_update_jacobian_3rd_derivatives( + CellSimilarity::none, + data_set, + data, + output_data.jacobian_3rd_derivatives); + maybe_update_jacobian_pushed_forward_3rd_derivatives( + CellSimilarity::none, + data_set, + data, + output_data.jacobian_pushed_forward_3rd_derivatives); + + maybe_compute_face_data(mapping, + cell, + face_no, + subface_no, + quadrature.size(), + quadrature.get_weights(), + data, + output_data); + } + + + + /** + * Implementation of MappingQGeneric::transform() for generic tensors. + */ + template + inline void + transform_fields( + const ArrayView> & input, + const MappingKind mapping_kind, + const typename Mapping::InternalDataBase &mapping_data, + const ArrayView> & output) + { + AssertDimension(input.size(), output.size()); + Assert((dynamic_cast::InternalData *>( + &mapping_data) != nullptr), + ExcInternalError()); + const typename dealii::MappingQGeneric::InternalData + &data = + static_cast:: + InternalData &>(mapping_data); + + switch (mapping_kind) + { + case mapping_contravariant: + { + Assert(data.update_each & update_contravariant_transformation, + typename FEValuesBase::ExcAccessToUninitializedField( + "update_contravariant_transformation")); + + for (unsigned int i = 0; i < output.size(); ++i) + output[i] = + apply_transformation(data.contravariant[i], input[i]); + + return; + } + + case mapping_piola: + { + Assert(data.update_each & update_contravariant_transformation, + typename FEValuesBase::ExcAccessToUninitializedField( + "update_contravariant_transformation")); + Assert(data.update_each & update_volume_elements, + typename FEValuesBase::ExcAccessToUninitializedField( + "update_volume_elements")); + Assert(rank == 1, ExcMessage("Only for rank 1")); + if (rank != 1) + return; + + for (unsigned int i = 0; i < output.size(); ++i) + { + output[i] = + apply_transformation(data.contravariant[i], input[i]); + output[i] /= data.volume_elements[i]; + } + return; + } + // We still allow this operation as in the + // reference cell Derivatives are Tensor + // rather than DerivativeForm + case mapping_covariant: + { + Assert(data.update_each & update_contravariant_transformation, + typename FEValuesBase::ExcAccessToUninitializedField( + "update_covariant_transformation")); + + for (unsigned int i = 0; i < output.size(); ++i) + output[i] = apply_transformation(data.covariant[i], input[i]); + + return; + } + + default: + Assert(false, ExcNotImplemented()); + } + } + + + + /** + * Implementation of MappingQGeneric::transform() for gradients. + */ + template + inline void + transform_gradients( + const ArrayView> & input, + const MappingKind mapping_kind, + const typename Mapping::InternalDataBase &mapping_data, + const ArrayView> & output) + { + AssertDimension(input.size(), output.size()); + Assert((dynamic_cast::InternalData *>( + &mapping_data) != nullptr), + ExcInternalError()); + const typename dealii::MappingQGeneric::InternalData + &data = + static_cast:: + InternalData &>(mapping_data); + + switch (mapping_kind) + { + case mapping_contravariant_gradient: + { + Assert(data.update_each & update_covariant_transformation, + typename FEValuesBase::ExcAccessToUninitializedField( + "update_covariant_transformation")); + Assert(data.update_each & update_contravariant_transformation, + typename FEValuesBase::ExcAccessToUninitializedField( + "update_contravariant_transformation")); + Assert(rank == 2, ExcMessage("Only for rank 2")); + + for (unsigned int i = 0; i < output.size(); ++i) + { + const DerivativeForm<1, spacedim, dim> A = + apply_transformation(data.contravariant[i], + transpose(input[i])); + output[i] = + apply_transformation(data.covariant[i], A.transpose()); + } + + return; + } + + case mapping_covariant_gradient: + { + Assert(data.update_each & update_covariant_transformation, + typename FEValuesBase::ExcAccessToUninitializedField( + "update_covariant_transformation")); + Assert(rank == 2, ExcMessage("Only for rank 2")); + + for (unsigned int i = 0; i < output.size(); ++i) + { + const DerivativeForm<1, spacedim, dim> A = + apply_transformation(data.covariant[i], + transpose(input[i])); + output[i] = + apply_transformation(data.covariant[i], A.transpose()); + } + + return; + } + + case mapping_piola_gradient: + { + Assert(data.update_each & update_covariant_transformation, + typename FEValuesBase::ExcAccessToUninitializedField( + "update_covariant_transformation")); + Assert(data.update_each & update_contravariant_transformation, + typename FEValuesBase::ExcAccessToUninitializedField( + "update_contravariant_transformation")); + Assert(data.update_each & update_volume_elements, + typename FEValuesBase::ExcAccessToUninitializedField( + "update_volume_elements")); + Assert(rank == 2, ExcMessage("Only for rank 2")); + + for (unsigned int i = 0; i < output.size(); ++i) + { + const DerivativeForm<1, spacedim, dim> A = + apply_transformation(data.covariant[i], input[i]); + const Tensor<2, spacedim> T = + apply_transformation(data.contravariant[i], A.transpose()); + + output[i] = transpose(T); + output[i] /= data.volume_elements[i]; + } + + return; + } + + default: + Assert(false, ExcNotImplemented()); + } + } + + + + /** + * Implementation of MappingQGeneric::transform() for hessians. + */ + template + inline void + transform_hessians( + const ArrayView> & input, + const MappingKind mapping_kind, + const typename Mapping::InternalDataBase &mapping_data, + const ArrayView> & output) + { + AssertDimension(input.size(), output.size()); + Assert((dynamic_cast::InternalData *>( + &mapping_data) != nullptr), + ExcInternalError()); + const typename dealii::MappingQGeneric::InternalData + &data = + static_cast:: + InternalData &>(mapping_data); + + switch (mapping_kind) + { + case mapping_contravariant_hessian: + { + Assert(data.update_each & update_covariant_transformation, + typename FEValuesBase::ExcAccessToUninitializedField( + "update_covariant_transformation")); + Assert(data.update_each & update_contravariant_transformation, + typename FEValuesBase::ExcAccessToUninitializedField( + "update_contravariant_transformation")); + + for (unsigned int q = 0; q < output.size(); ++q) + for (unsigned int i = 0; i < spacedim; ++i) + { + double tmp1[dim][dim]; + for (unsigned int J = 0; J < dim; ++J) + for (unsigned int K = 0; K < dim; ++K) + { + tmp1[J][K] = + data.contravariant[q][i][0] * input[q][0][J][K]; + for (unsigned int I = 1; I < dim; ++I) + tmp1[J][K] += + data.contravariant[q][i][I] * input[q][I][J][K]; + } + for (unsigned int j = 0; j < spacedim; ++j) + { + double tmp2[dim]; + for (unsigned int K = 0; K < dim; ++K) + { + tmp2[K] = data.covariant[q][j][0] * tmp1[0][K]; + for (unsigned int J = 1; J < dim; ++J) + tmp2[K] += data.covariant[q][j][J] * tmp1[J][K]; + } + for (unsigned int k = 0; k < spacedim; ++k) + { + output[q][i][j][k] = + data.covariant[q][k][0] * tmp2[0]; + for (unsigned int K = 1; K < dim; ++K) + output[q][i][j][k] += + data.covariant[q][k][K] * tmp2[K]; + } + } + } + return; + } + + case mapping_covariant_hessian: + { + Assert(data.update_each & update_covariant_transformation, + typename FEValuesBase::ExcAccessToUninitializedField( + "update_covariant_transformation")); + + for (unsigned int q = 0; q < output.size(); ++q) + for (unsigned int i = 0; i < spacedim; ++i) + { + double tmp1[dim][dim]; + for (unsigned int J = 0; J < dim; ++J) + for (unsigned int K = 0; K < dim; ++K) + { + tmp1[J][K] = + data.covariant[q][i][0] * input[q][0][J][K]; + for (unsigned int I = 1; I < dim; ++I) + tmp1[J][K] += + data.covariant[q][i][I] * input[q][I][J][K]; + } + for (unsigned int j = 0; j < spacedim; ++j) + { + double tmp2[dim]; + for (unsigned int K = 0; K < dim; ++K) + { + tmp2[K] = data.covariant[q][j][0] * tmp1[0][K]; + for (unsigned int J = 1; J < dim; ++J) + tmp2[K] += data.covariant[q][j][J] * tmp1[J][K]; + } + for (unsigned int k = 0; k < spacedim; ++k) + { + output[q][i][j][k] = + data.covariant[q][k][0] * tmp2[0]; + for (unsigned int K = 1; K < dim; ++K) + output[q][i][j][k] += + data.covariant[q][k][K] * tmp2[K]; + } + } + } + + return; + } + + case mapping_piola_hessian: + { + Assert(data.update_each & update_covariant_transformation, + typename FEValuesBase::ExcAccessToUninitializedField( + "update_covariant_transformation")); + Assert(data.update_each & update_contravariant_transformation, + typename FEValuesBase::ExcAccessToUninitializedField( + "update_contravariant_transformation")); + Assert(data.update_each & update_volume_elements, + typename FEValuesBase::ExcAccessToUninitializedField( + "update_volume_elements")); + + for (unsigned int q = 0; q < output.size(); ++q) + for (unsigned int i = 0; i < spacedim; ++i) + { + double factor[dim]; + for (unsigned int I = 0; I < dim; ++I) + factor[I] = + data.contravariant[q][i][I] / data.volume_elements[q]; + double tmp1[dim][dim]; + for (unsigned int J = 0; J < dim; ++J) + for (unsigned int K = 0; K < dim; ++K) + { + tmp1[J][K] = factor[0] * input[q][0][J][K]; + for (unsigned int I = 1; I < dim; ++I) + tmp1[J][K] += factor[I] * input[q][I][J][K]; + } + for (unsigned int j = 0; j < spacedim; ++j) + { + double tmp2[dim]; + for (unsigned int K = 0; K < dim; ++K) + { + tmp2[K] = data.covariant[q][j][0] * tmp1[0][K]; + for (unsigned int J = 1; J < dim; ++J) + tmp2[K] += data.covariant[q][j][J] * tmp1[J][K]; + } + for (unsigned int k = 0; k < spacedim; ++k) + { + output[q][i][j][k] = + data.covariant[q][k][0] * tmp2[0]; + for (unsigned int K = 1; K < dim; ++K) + output[q][i][j][k] += + data.covariant[q][k][K] * tmp2[K]; + } + } + } + + return; + } + + default: + Assert(false, ExcNotImplemented()); + } + } + + + + /** + * Implementation of MappingQGeneric::transform() for DerivativeForm + * arguments. + */ + template + inline void + transform_differential_forms( + const ArrayView> &input, + const MappingKind mapping_kind, + const typename Mapping::InternalDataBase & mapping_data, + const ArrayView> & output) + { + AssertDimension(input.size(), output.size()); + Assert((dynamic_cast::InternalData *>( + &mapping_data) != nullptr), + ExcInternalError()); + const typename dealii::MappingQGeneric::InternalData + &data = + static_cast:: + InternalData &>(mapping_data); + + switch (mapping_kind) + { + case mapping_covariant: + { + Assert(data.update_each & update_contravariant_transformation, + typename FEValuesBase::ExcAccessToUninitializedField( + "update_covariant_transformation")); + + for (unsigned int i = 0; i < output.size(); ++i) + output[i] = apply_transformation(data.covariant[i], input[i]); + + return; + } + default: + Assert(false, ExcNotImplemented()); + } + } + } // namespace MappingQGenericImplementation +} // namespace internal + +DEAL_II_NAMESPACE_CLOSE + +#endif diff --git a/source/fe/mapping_q_generic.cc b/source/fe/mapping_q_generic.cc index 4e8b7f0cf7..f452d68872 100644 --- a/source/fe/mapping_q_generic.cc +++ b/source/fe/mapping_q_generic.cc @@ -28,6 +28,7 @@ #include #include #include +#include #include #include @@ -35,11 +36,6 @@ #include #include -#include - -#include -#include -#include #include @@ -53,161 +49,6 @@ 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::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::array, GeometryInfo<2>::vertices_per_cell> - & vertices, - const Point &p) - { - Assert(spacedim == 2, ExcInternalError()); - - // For accuracy reasons, we do all arithmetic in extended precision - // (long double). This has a noticeable effect on the hit rate for - // borderline cases and thus makes the algorithm more robust. - const long double x = p(0); - const long double y = p(1); - - const long double x0 = vertices[0](0); - const long double x1 = vertices[1](0); - const long double x2 = vertices[2](0); - const long double x3 = vertices[3](0); - - const long double y0 = vertices[0](1); - const long double y1 = vertices[1](1); - const long double y2 = vertices[2](1); - const long double y3 = vertices[3](1); - - const long double a = (x1 - x3) * (y0 - y2) - (x0 - x2) * (y1 - y3); - const long double b = -(x0 - x1 - x2 + x3) * y + - (x - 2 * x1 + x3) * y0 - (x - 2 * x0 + x2) * y1 - - (x - x1) * y2 + (x - x0) * y3; - const long double c = (x0 - x1) * y - (x - x1) * y0 + (x - x0) * y1; - - const long 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) - AssertThrow( - discriminant > 0.0, - (typename Mapping::ExcTransformationFailed())); - - long double eta1; - long double eta2; - const long double sqrt_discriminant = std::sqrt(discriminant); - // special case #1: if a is near-zero to make the discriminant exactly - // equal b, then use the linear formula - if (b != 0.0 && std::abs(b) == sqrt_discriminant) - { - eta1 = -c / b; - eta2 = -c / b; - } - // special case #2: a is zero for parallelograms and very small for - // near-parallelograms: - else if (std::abs(a) < 1e-8 * std::abs(b)) - { - // if both a and c are very small then the root should be near - // zero: this first case will capture that - eta1 = 2 * c / (-b - sqrt_discriminant); - eta2 = 2 * c / (-b + sqrt_discriminant); - } - // finally, use the plain version: - else - { - eta1 = (-b - sqrt_discriminant) / (2 * a); - eta2 = (-b + sqrt_discriminant) / (2 * a); - } - // pick the one closer to the center of the cell. - const long 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 long double subexpr0 = -eta * x2 + x0 * (eta - 1); - const long double xi_denominator0 = - eta * x3 - x1 * (eta - 1) + subexpr0; - const long 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 {xi, static_cast(eta)}; - } - else - { - const long double max_y = - std::max(std::max(std::abs(y0), std::abs(y1)), - std::max(std::abs(y2), std::abs(y3))); - const long double subexpr1 = -eta * y2 + y0 * (eta - 1); - const long 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 {xi, static_cast(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 {std::numeric_limits::quiet_NaN(), - std::numeric_limits::quiet_NaN()}; - } - - - - template - Point<3> - transform_real_to_unit_cell( - const std::array, GeometryInfo<3>::vertices_per_cell> - & /*vertices*/, - const Point & /*p*/) - { - // It should not be possible to get here - Assert(false, ExcInternalError()); - return Point<3>(); - } - } // namespace - } // namespace MappingQ1 -} // namespace internal - - - template MappingQGeneric::InternalData::InternalData( const unsigned int polynomial_degree) @@ -239,6 +80,7 @@ MappingQGeneric::InternalData::memory_consumption() const } + template void MappingQGeneric::InternalData::initialize( @@ -517,1387 +359,6 @@ MappingQGeneric::InternalData::compute_shape_function_values( -namespace internal -{ - namespace MappingQGenericImplementation - { - namespace - { - /** - * This function is needed by the constructor of - * MappingQ for dim= 2 and 3. - * - * For the definition of the @p support_point_weights_on_quad please - * refer to the description of TransfiniteInterpolationManifold. - */ - dealii::Table<2, double> - compute_support_point_weights_on_quad( - const unsigned int polynomial_degree) - { - dealii::Table<2, double> loqvs; - - // we are asked to compute weights for interior support points, but - // there are no interior points if degree==1 - if (polynomial_degree == 1) - return loqvs; - - const unsigned int M = polynomial_degree - 1; - const unsigned int n_inner_2d = M * M; - const unsigned int n_outer_2d = 4 + 4 * M; - - // set the weights of transfinite interpolation - loqvs.reinit(n_inner_2d, n_outer_2d); - QGaussLobatto<2> gl(polynomial_degree + 1); - for (unsigned int i = 0; i < M; ++i) - for (unsigned int j = 0; j < M; ++j) - { - const Point<2> p = - gl.point((i + 1) * (polynomial_degree + 1) + (j + 1)); - const unsigned int index_table = i * M + j; - for (unsigned int v = 0; v < 4; ++v) - loqvs(index_table, v) = - -GeometryInfo<2>::d_linear_shape_function(p, v); - loqvs(index_table, 4 + i) = 1. - p[0]; - loqvs(index_table, 4 + i + M) = p[0]; - loqvs(index_table, 4 + j + 2 * M) = 1. - p[1]; - loqvs(index_table, 4 + j + 3 * M) = p[1]; - } - - // the sum of weights of the points at the outer rim should be one. - // check this - for (unsigned int unit_point = 0; unit_point < n_inner_2d; ++unit_point) - Assert(std::fabs(std::accumulate(loqvs[unit_point].begin(), - loqvs[unit_point].end(), - 0.) - - 1) < 1e-13 * polynomial_degree, - ExcInternalError()); - - return loqvs; - } - - - - /** - * This function is needed by the constructor of MappingQ<3>. - * - * For the definition of the @p support_point_weights_on_quad please - * refer to the description of TransfiniteInterpolationManifold. - */ - dealii::Table<2, double> - compute_support_point_weights_on_hex(const unsigned int polynomial_degree) - { - dealii::Table<2, double> lohvs; - - // we are asked to compute weights for interior support points, but - // there are no interior points if degree==1 - if (polynomial_degree == 1) - return lohvs; - - const unsigned int M = polynomial_degree - 1; - - const unsigned int n_inner = Utilities::fixed_power<3>(M); - const unsigned int n_outer = 8 + 12 * M + 6 * M * M; - - // set the weights of transfinite interpolation - lohvs.reinit(n_inner, n_outer); - QGaussLobatto<3> gl(polynomial_degree + 1); - for (unsigned int i = 0; i < M; ++i) - for (unsigned int j = 0; j < M; ++j) - for (unsigned int k = 0; k < M; ++k) - { - const Point<3> p = gl.point((i + 1) * (M + 2) * (M + 2) + - (j + 1) * (M + 2) + (k + 1)); - const unsigned int index_table = i * M * M + j * M + k; - - // vertices - for (unsigned int v = 0; v < 8; ++v) - lohvs(index_table, v) = - GeometryInfo<3>::d_linear_shape_function(p, v); - - // lines - { - constexpr std::array line_coordinates_y( - {{0, 1, 4, 5}}); - const Point<2> py(p[0], p[2]); - for (unsigned int l = 0; l < 4; ++l) - lohvs(index_table, 8 + line_coordinates_y[l] * M + j) = - -GeometryInfo<2>::d_linear_shape_function(py, l); - } - - { - constexpr std::array line_coordinates_x( - {{2, 3, 6, 7}}); - const Point<2> px(p[1], p[2]); - for (unsigned int l = 0; l < 4; ++l) - lohvs(index_table, 8 + line_coordinates_x[l] * M + k) = - -GeometryInfo<2>::d_linear_shape_function(px, l); - } - - { - constexpr std::array line_coordinates_z( - {{8, 9, 10, 11}}); - const Point<2> pz(p[0], p[1]); - for (unsigned int l = 0; l < 4; ++l) - lohvs(index_table, 8 + line_coordinates_z[l] * M + i) = - -GeometryInfo<2>::d_linear_shape_function(pz, l); - } - - // quads - lohvs(index_table, 8 + 12 * M + 0 * M * M + i * M + j) = - 1. - p[0]; - lohvs(index_table, 8 + 12 * M + 1 * M * M + i * M + j) = p[0]; - lohvs(index_table, 8 + 12 * M + 2 * M * M + k * M + i) = - 1. - p[1]; - lohvs(index_table, 8 + 12 * M + 3 * M * M + k * M + i) = p[1]; - lohvs(index_table, 8 + 12 * M + 4 * M * M + j * M + k) = - 1. - p[2]; - lohvs(index_table, 8 + 12 * M + 5 * M * M + j * M + k) = p[2]; - } - - // the sum of weights of the points at the outer rim should be one. - // check this - for (unsigned int unit_point = 0; unit_point < n_inner; ++unit_point) - Assert(std::fabs(std::accumulate(lohvs[unit_point].begin(), - lohvs[unit_point].end(), - 0.) - - 1) < 1e-13 * polynomial_degree, - ExcInternalError()); - - return lohvs; - } - - - - /** - * This function collects the output of - * compute_support_point_weights_on_{quad,hex} in a single data structure. - */ - std::vector> - compute_support_point_weights_perimeter_to_interior( - const unsigned int polynomial_degree, - const unsigned int dim) - { - Assert(dim > 0 && dim <= 3, ExcImpossibleInDim(dim)); - std::vector> output(dim); - if (polynomial_degree <= 1) - return output; - - // fill the 1D interior weights - QGaussLobatto<1> quadrature(polynomial_degree + 1); - output[0].reinit(polynomial_degree - 1, - GeometryInfo<1>::vertices_per_cell); - for (unsigned int q = 0; q < polynomial_degree - 1; ++q) - for (const unsigned int i : GeometryInfo<1>::vertex_indices()) - output[0](q, i) = - GeometryInfo<1>::d_linear_shape_function(quadrature.point(q + 1), - i); - - if (dim > 1) - output[1] = compute_support_point_weights_on_quad(polynomial_degree); - - if (dim > 2) - output[2] = compute_support_point_weights_on_hex(polynomial_degree); - - return output; - } - - /** - * Collects all interior points for the various dimensions. - */ - template - dealii::Table<2, double> - compute_support_point_weights_cell(const unsigned int polynomial_degree) - { - Assert(dim > 0 && dim <= 3, ExcImpossibleInDim(dim)); - if (polynomial_degree <= 1) - return dealii::Table<2, double>(); - - QGaussLobatto quadrature(polynomial_degree + 1); - const std::vector h2l = - FETools::hierarchic_to_lexicographic_numbering( - polynomial_degree); - - dealii::Table<2, double> output(quadrature.size() - - GeometryInfo::vertices_per_cell, - GeometryInfo::vertices_per_cell); - for (unsigned int q = 0; q < output.size(0); ++q) - for (const unsigned int i : GeometryInfo::vertex_indices()) - output(q, i) = GeometryInfo::d_linear_shape_function( - quadrature.point(h2l[q + GeometryInfo::vertices_per_cell]), - i); - - return output; - } - - - - /** - * Using the relative weights of the shape functions evaluated at - * one point on the reference cell (and stored in data.shape_values - * and accessed via data.shape(0,i)) and the locations of mapping - * support points (stored in data.mapping_support_points), compute - * the mapped location of that point in real space. - */ - template - Point - compute_mapped_location_of_point( - const typename dealii::MappingQGeneric::InternalData - &data) - { - AssertDimension(data.shape_values.size(), - data.mapping_support_points.size()); - - // use now the InternalData to compute the point in real space. - Point p_real; - for (unsigned int i = 0; i < data.mapping_support_points.size(); ++i) - p_real += data.mapping_support_points[i] * data.shape(0, i); - - return p_real; - } - - - /** - * Implementation of transform_real_to_unit_cell for either type double - * or VectorizedArray - */ - template - Point - do_transform_real_to_unit_cell_internal( - const Point & p, - const Point & initial_p_unit, - const std::vector> & points, - const std::vector> &polynomials_1d, - const std::vector & renumber) - { - AssertDimension(points.size(), - Utilities::pow(polynomials_1d.size(), dim)); - - // Newton iteration to solve - // f(x)=p(x)-p=0 - // where we are looking for 'x' and p(x) is the forward transformation - // from unit to real cell. We solve this using a Newton iteration - // x_{n+1}=x_n-[f'(x)]^{-1}f(x) - // The start value is set to be the linear approximation to the cell - - // The shape values and derivatives of the mapping at this point are - // previously computed. - - Point p_unit = initial_p_unit; - auto p_real = internal::evaluate_tensor_product_value_and_gradient( - polynomials_1d, points, p_unit, polynomials_1d.size() == 2, renumber); - - Tensor<1, spacedim, Number> f = p_real.first - p; - - // early out if we already have our point in all SIMD lanes, i.e., - // f.norm_square() < 1e-24 * p_real.second[0].norm_square(). To enable - // this step also for VectorizedArray where we do not have operator <, - // compare the result to zero which is defined for SIMD types - if (std::max(Number(0.), - f.norm_square() - - 1e-24 * p_real.second[0].norm_square()) == Number(0.)) - return p_unit; - - // we need to compare the position of the computed p(x) against the - // given point 'p'. We will terminate the iteration and return 'x' if - // they are less than eps apart. The question is how to choose eps -- - // or, put maybe more generally: in which norm we want these 'p' and - // 'p(x)' to be eps apart. - // - // the question is difficult since we may have to deal with very - // elongated cells where we may achieve 1e-12*h for the distance of - // these two points in the 'long' direction, but achieving this - // tolerance in the 'short' direction of the cell may not be possible - // - // what we do instead is then to terminate iterations if - // \| p(x) - p \|_A < eps - // where the A-norm is somehow induced by the transformation of the - // cell. in particular, we want to measure distances relative to the - // sizes of the cell in its principal directions. - // - // to define what exactly A should be, note that to first order we have - // the following (assuming that x* is the solution of the problem, i.e., - // p(x*)=p): - // p(x) - p = p(x) - p(x*) - // = -grad p(x) * (x*-x) + higher order terms - // This suggest to measure with a norm that corresponds to - // A = {[grad p(x]^T [grad p(x)]}^{-1} - // because then - // \| p(x) - p \|_A \approx \| x - x* \| - // Consequently, we will try to enforce that - // \| p(x) - p \|_A = \| f \| <= eps - // - // Note that using this norm is a bit dangerous since the norm changes - // in every iteration (A isn't fixed by depending on xk). However, if - // the cell is not too deformed (it may be stretched, but not twisted) - // then the mapping is almost linear and A is indeed constant or - // nearly so. - const double eps = 1.e-11; - const unsigned int newton_iteration_limit = 20; - - Point invalid_point; - invalid_point[0] = std::numeric_limits::infinity(); - bool try_project_to_unit_cell = false; - - unsigned int newton_iteration = 0; - Number f_weighted_norm_square = 1.; - Number last_f_weighted_norm_square = 1.; - - do - { -#ifdef DEBUG_TRANSFORM_REAL_TO_UNIT_CELL - std::cout << "Newton iteration " << newton_iteration - << " point guess " << p_unit << std::endl; -#endif - - // f'(x) - Tensor<2, spacedim, Number> df; - for (unsigned int d = 0; d < spacedim; ++d) - for (unsigned int e = 0; e < dim; ++e) - df[d][e] = p_real.second[e][d]; - - // check determinand(df) > 0 on all SIMD lanes - if (!(std::min(determinant(df), - Number(std::numeric_limits::min())) == - Number(std::numeric_limits::min()))) - { - // We allow to enter this function with an initial guess - // outside the unit cell. We might have run into invalid - // Jacobians due to that, so we should at least try once to go - // back to the unit cell and go on with the Newton iteration - // from there. Since the outside case is unlikely, we can - // afford spending the extra effort at this place. - if (try_project_to_unit_cell == false) - { - p_unit = GeometryInfo::project_to_unit_cell(p_unit); - p_real = - internal::evaluate_tensor_product_value_and_gradient( - polynomials_1d, - points, - p_unit, - polynomials_1d.size() == 2, - renumber); - f = p_real.first - p; - f_weighted_norm_square = 1.; - last_f_weighted_norm_square = 1; - try_project_to_unit_cell = true; - continue; - } - else - return invalid_point; - } - - // Solve [f'(x)]d=f(x) - const Tensor<2, spacedim, Number> df_inverse = invert(df); - const Tensor<1, spacedim, Number> delta = df_inverse * f; - last_f_weighted_norm_square = delta.norm_square(); - -#ifdef DEBUG_TRANSFORM_REAL_TO_UNIT_CELL - std::cout << " delta=" << delta << std::endl; -#endif - - // do a line search - double step_length = 1; - do - { - // update of p_unit. The spacedim-th component of transformed - // point is simply ignored in codimension one case. When this - // component is not zero, then we are projecting the point to - // the surface or curve identified by the cell. - Point p_unit_trial = p_unit; - for (unsigned int i = 0; i < dim; ++i) - p_unit_trial[i] -= step_length * delta[i]; - - // shape values and derivatives at new p_unit point - const auto p_real_trial = - internal::evaluate_tensor_product_value_and_gradient( - polynomials_1d, - points, - p_unit_trial, - polynomials_1d.size() == 2, - renumber); - const Tensor<1, spacedim, Number> f_trial = - p_real_trial.first - p; - f_weighted_norm_square = (df_inverse * f_trial).norm_square(); - -#ifdef DEBUG_TRANSFORM_REAL_TO_UNIT_CELL - std::cout << " step_length=" << step_length << std::endl - << " ||f || =" << f.norm() << std::endl - << " ||f*|| =" << f_trial.norm() << std::endl - << " ||f*||_A =" - << std::sqrt(f_weighted_norm_square) << std::endl; -#endif - - // See if we are making progress with the current step length - // and if not, reduce it by a factor of two and try again. - // - // Strictly speaking, we should probably use the same norm as we - // use for the outer algorithm. In practice, line search is just - // a crutch to find a "reasonable" step length, and so using the - // l2 norm is probably just fine. - // - // check f_trial.norm() < f.norm() in SIMD form. This is a bit - // more complicated because some SIMD lanes might not be doing - // any progress any more as they have already reached roundoff - // accuracy. We define that as the case of updates less than - // 1e-6. The tolerance might seem coarse but since we are - // dealing with a Newton iteration of a polynomial function we - // either converge quadratically or not any more. Thus, our - // selection is to terminate if either the norm of f is - // decreasing or that threshold of 1e-6 is reached. - if (std::max(f_weighted_norm_square - 1e-6 * 1e-6, Number(0.)) * - std::max(f_trial.norm_square() - f.norm_square(), - Number(0.)) == - Number(0.)) - { - p_real = p_real_trial; - p_unit = p_unit_trial; - f = f_trial; - break; - } - else if (step_length > 0.05) - step_length *= 0.5; - else - break; - } - while (true); - - // In case we terminated the line search due to the step becoming - // too small, we give the iteration another try with the - // projection of the initial guess to the unit cell before we give - // up, just like for the negative determinant case. - if (step_length <= 0.05 && try_project_to_unit_cell == false) - { - p_unit = GeometryInfo::project_to_unit_cell(p_unit); - p_real = internal::evaluate_tensor_product_value_and_gradient( - polynomials_1d, - points, - p_unit, - polynomials_1d.size() == 2, - renumber); - f = p_real.first - p; - f_weighted_norm_square = 1.; - last_f_weighted_norm_square = 1; - try_project_to_unit_cell = true; - continue; - } - else if (step_length <= 0.05) - return invalid_point; - - ++newton_iteration; - if (newton_iteration > newton_iteration_limit) - return invalid_point; - } - // Stop if f_weighted_norm_square <= eps^2 on all SIMD lanes or if the - // weighted norm is less than 1e-6 and the improvement against the - // previous step was less than a factor of 10 (in that regime, we - // either have quadratic convergence or roundoff errors due to a bad - // mapping) - while ( - !(std::max(f_weighted_norm_square - eps * eps, Number(0.)) * - std::max(last_f_weighted_norm_square - - std::min(f_weighted_norm_square, Number(1e-6 * 1e-6)) * - 100., - Number(0.)) == - Number(0.))); - - 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 dealii::Triangulation::cell_iterator &cell, - const Point & p, - const Point &initial_p_unit, - typename dealii::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 < 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 &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(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 &>(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>(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 &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(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::ExcTransformationFailed())); - - return p_unit; - } - - /** - * In case the quadrature formula is a tensor product, this is a - * replacement for maybe_compute_q_points(), maybe_update_Jacobians() and - * maybe_update_jacobian_grads() - */ - template - void - maybe_update_q_points_Jacobians_and_grads_tensor( - const CellSimilarity::Similarity cell_similarity, - const typename dealii::MappingQGeneric::InternalData - & data, - std::vector> & quadrature_points, - std::vector> &jacobian_grads) - { - const UpdateFlags update_flags = data.update_each; - - const unsigned int n_shape_values = data.n_shape_functions; - const unsigned int n_q_points = data.shape_info.n_q_points; - constexpr unsigned int n_lanes = VectorizedArray::size(); - constexpr unsigned int n_comp = 1 + (spacedim - 1) / n_lanes; - constexpr unsigned int n_hessians = (dim * (dim + 1)) / 2; - - EvaluationFlags::EvaluationFlags evaluation_flag = - (update_flags & update_quadrature_points ? EvaluationFlags::values : - EvaluationFlags::nothing) | - ((cell_similarity != CellSimilarity::translation) && - (update_flags & update_contravariant_transformation) ? - EvaluationFlags::gradients : - EvaluationFlags::nothing) | - ((cell_similarity != CellSimilarity::translation) && - (update_flags & update_jacobian_grads) ? - EvaluationFlags::hessians : - EvaluationFlags::nothing); - - Assert(!(evaluation_flag & EvaluationFlags::values) || n_q_points > 0, - ExcInternalError()); - Assert(!(evaluation_flag & EvaluationFlags::values) || - n_q_points == quadrature_points.size(), - ExcDimensionMismatch(n_q_points, quadrature_points.size())); - Assert(!(evaluation_flag & EvaluationFlags::gradients) || - data.n_shape_functions > 0, - ExcInternalError()); - Assert(!(evaluation_flag & EvaluationFlags::gradients) || - n_q_points == data.contravariant.size(), - ExcDimensionMismatch(n_q_points, data.contravariant.size())); - Assert(!(evaluation_flag & EvaluationFlags::hessians) || - n_q_points == jacobian_grads.size(), - ExcDimensionMismatch(n_q_points, jacobian_grads.size())); - - // shortcut in case we have an identity interpolation and only request - // the quadrature points - if (evaluation_flag == EvaluationFlags::values && - data.shape_info.element_type == - internal::MatrixFreeFunctions::tensor_symmetric_collocation) - { - for (unsigned int q = 0; q < n_q_points; ++q) - quadrature_points[q] = - data.mapping_support_points[data.shape_info - .lexicographic_numbering[q]]; - return; - } - - // prepare arrays - if (evaluation_flag != EvaluationFlags::nothing) - { - data.values_dofs.resize(n_comp * n_shape_values); - data.values_quad.resize(n_comp * n_q_points); - data.gradients_quad.resize(n_comp * n_q_points * dim); - data.scratch.resize(2 * std::max(n_q_points, n_shape_values)); - - if (evaluation_flag & EvaluationFlags::hessians) - data.hessians_quad.resize(n_comp * n_q_points * n_hessians); - - const std::vector &renumber_to_lexicographic = - data.shape_info.lexicographic_numbering; - for (unsigned int i = 0; i < n_shape_values; ++i) - for (unsigned int d = 0; d < spacedim; ++d) - { - const unsigned int in_comp = d % n_lanes; - const unsigned int out_comp = d / n_lanes; - data.values_dofs[out_comp * n_shape_values + i][in_comp] = - data - .mapping_support_points[renumber_to_lexicographic[i]][d]; - } - - // do the actual tensorized evaluation - internal::FEEvaluationFactory< - dim, - double, - VectorizedArray>::evaluate(n_comp, - evaluation_flag, - data.shape_info, - data.values_dofs.begin(), - data.values_quad.begin(), - data.gradients_quad.begin(), - data.hessians_quad.begin(), - data.scratch.begin()); - } - - // do the postprocessing - if (evaluation_flag & EvaluationFlags::values) - { - for (unsigned int out_comp = 0; out_comp < n_comp; ++out_comp) - for (unsigned int i = 0; i < n_q_points; ++i) - for (unsigned int in_comp = 0; - in_comp < n_lanes && - in_comp < spacedim - out_comp * n_lanes; - ++in_comp) - quadrature_points[i][out_comp * n_lanes + in_comp] = - data.values_quad[out_comp * n_q_points + i][in_comp]; - } - - if (evaluation_flag & EvaluationFlags::gradients) - { - std::fill(data.contravariant.begin(), - data.contravariant.end(), - DerivativeForm<1, dim, spacedim>()); - // We need to reinterpret the data after evaluate has been applied. - for (unsigned int out_comp = 0; out_comp < n_comp; ++out_comp) - for (unsigned int point = 0; point < n_q_points; ++point) - for (unsigned int j = 0; j < dim; ++j) - for (unsigned int in_comp = 0; - in_comp < n_lanes && - in_comp < spacedim - out_comp * n_lanes; - ++in_comp) - { - const unsigned int total_number = point * dim + j; - const unsigned int new_comp = total_number / n_q_points; - const unsigned int new_point = total_number % n_q_points; - data.contravariant[new_point][out_comp * n_lanes + - in_comp][new_comp] = - data.gradients_quad[(out_comp * n_q_points + point) * - dim + - j][in_comp]; - } - } - if (update_flags & update_covariant_transformation) - if (cell_similarity != CellSimilarity::translation) - for (unsigned int point = 0; point < n_q_points; ++point) - data.covariant[point] = - (data.contravariant[point]).covariant_form(); - - if (update_flags & update_volume_elements) - if (cell_similarity != CellSimilarity::translation) - for (unsigned int point = 0; point < n_q_points; ++point) - data.volume_elements[point] = - data.contravariant[point].determinant(); - - if (evaluation_flag & EvaluationFlags::hessians) - { - constexpr int desymmetrize_3d[6][2] = { - {0, 0}, {1, 1}, {2, 2}, {0, 1}, {0, 2}, {1, 2}}; - constexpr int desymmetrize_2d[3][2] = {{0, 0}, {1, 1}, {0, 1}}; - - // We need to reinterpret the data after evaluate has been applied. - for (unsigned int out_comp = 0; out_comp < n_comp; ++out_comp) - for (unsigned int point = 0; point < n_q_points; ++point) - for (unsigned int j = 0; j < n_hessians; ++j) - for (unsigned int in_comp = 0; - in_comp < n_lanes && - in_comp < spacedim - out_comp * n_lanes; - ++in_comp) - { - const unsigned int total_number = point * n_hessians + j; - const unsigned int new_point = total_number % n_q_points; - const unsigned int new_hessian_comp = - total_number / n_q_points; - const unsigned int new_hessian_comp_i = - dim == 2 ? desymmetrize_2d[new_hessian_comp][0] : - desymmetrize_3d[new_hessian_comp][0]; - const unsigned int new_hessian_comp_j = - dim == 2 ? desymmetrize_2d[new_hessian_comp][1] : - desymmetrize_3d[new_hessian_comp][1]; - const double value = - data.hessians_quad[(out_comp * n_q_points + point) * - n_hessians + - j][in_comp]; - jacobian_grads[new_point][out_comp * n_lanes + in_comp] - [new_hessian_comp_i][new_hessian_comp_j] = - value; - jacobian_grads[new_point][out_comp * n_lanes + in_comp] - [new_hessian_comp_j][new_hessian_comp_i] = - value; - } - } - } - - - /** - * Compute the locations of quadrature points on the object described by - * the first argument (and the cell for which the mapping support points - * have already been set), but only if the update_flags of the @p data - * argument indicate so. - */ - template - void - maybe_compute_q_points( - const typename QProjector::DataSetDescriptor data_set, - const typename dealii::MappingQGeneric::InternalData - & data, - std::vector> &quadrature_points) - { - const UpdateFlags update_flags = data.update_each; - - if (update_flags & update_quadrature_points) - for (unsigned int point = 0; point < quadrature_points.size(); - ++point) - { - const double * shape = &data.shape(point + data_set, 0); - Point result = - (shape[0] * data.mapping_support_points[0]); - for (unsigned int k = 1; k < data.n_shape_functions; ++k) - for (unsigned int i = 0; i < spacedim; ++i) - result[i] += shape[k] * data.mapping_support_points[k][i]; - quadrature_points[point] = result; - } - } - - - - /** - * Update the co- and contravariant matrices as well as their determinant, - * for the cell - * described stored in the data object, but only if the update_flags of the @p data - * argument indicate so. - * - * Skip the computation if possible as indicated by the first argument. - */ - template - void - maybe_update_Jacobians( - const CellSimilarity::Similarity cell_similarity, - const typename dealii::QProjector::DataSetDescriptor data_set, - const typename dealii::MappingQGeneric::InternalData - &data) - { - const UpdateFlags update_flags = data.update_each; - - if (update_flags & update_contravariant_transformation) - // if the current cell is just a - // translation of the previous one, no - // need to recompute jacobians... - if (cell_similarity != CellSimilarity::translation) - { - const unsigned int n_q_points = data.contravariant.size(); - - std::fill(data.contravariant.begin(), - data.contravariant.end(), - DerivativeForm<1, dim, spacedim>()); - - Assert(data.n_shape_functions > 0, ExcInternalError()); - - const Tensor<1, spacedim> *supp_pts = - data.mapping_support_points.data(); - - for (unsigned int point = 0; point < n_q_points; ++point) - { - const Tensor<1, dim> *data_derv = - &data.derivative(point + data_set, 0); - - double result[spacedim][dim]; - - // peel away part of sum to avoid zeroing the - // entries and adding for the first time - for (unsigned int i = 0; i < spacedim; ++i) - for (unsigned int j = 0; j < dim; ++j) - result[i][j] = data_derv[0][j] * supp_pts[0][i]; - for (unsigned int k = 1; k < data.n_shape_functions; ++k) - for (unsigned int i = 0; i < spacedim; ++i) - for (unsigned int j = 0; j < dim; ++j) - result[i][j] += data_derv[k][j] * supp_pts[k][i]; - - // write result into contravariant data. for - // j=dim in the case dim - void - maybe_update_jacobian_grads( - const CellSimilarity::Similarity cell_similarity, - const typename QProjector::DataSetDescriptor data_set, - const typename dealii::MappingQGeneric::InternalData - & data, - std::vector> &jacobian_grads) - { - const UpdateFlags update_flags = data.update_each; - if (update_flags & update_jacobian_grads) - { - const unsigned int n_q_points = jacobian_grads.size(); - - if (cell_similarity != CellSimilarity::translation) - for (unsigned int point = 0; point < n_q_points; ++point) - { - const Tensor<2, dim> *second = - &data.second_derivative(point + data_set, 0); - double result[spacedim][dim][dim]; - for (unsigned int i = 0; i < spacedim; ++i) - for (unsigned int j = 0; j < dim; ++j) - for (unsigned int l = 0; l < dim; ++l) - result[i][j][l] = - (second[0][j][l] * data.mapping_support_points[0][i]); - for (unsigned int k = 1; k < data.n_shape_functions; ++k) - for (unsigned int i = 0; i < spacedim; ++i) - for (unsigned int j = 0; j < dim; ++j) - for (unsigned int l = 0; l < dim; ++l) - result[i][j][l] += - (second[k][j][l] * - data.mapping_support_points[k][i]); - - for (unsigned int i = 0; i < spacedim; ++i) - for (unsigned int j = 0; j < dim; ++j) - for (unsigned int l = 0; l < dim; ++l) - jacobian_grads[point][i][j][l] = result[i][j][l]; - } - } - } - - /** - * Update the Hessian of the transformation from unit to real cell, the - * Jacobian gradients, pushed forward to the real cell coordinates. - * - * Skip the computation if possible as indicated by the first argument. - */ - template - void - maybe_update_jacobian_pushed_forward_grads( - const CellSimilarity::Similarity cell_similarity, - const typename QProjector::DataSetDescriptor data_set, - const typename dealii::MappingQGeneric::InternalData - & data, - std::vector> &jacobian_pushed_forward_grads) - { - const UpdateFlags update_flags = data.update_each; - if (update_flags & update_jacobian_pushed_forward_grads) - { - const unsigned int n_q_points = - jacobian_pushed_forward_grads.size(); - - if (cell_similarity != CellSimilarity::translation) - { - double tmp[spacedim][spacedim][spacedim]; - for (unsigned int point = 0; point < n_q_points; ++point) - { - const Tensor<2, dim> *second = - &data.second_derivative(point + data_set, 0); - double result[spacedim][dim][dim]; - for (unsigned int i = 0; i < spacedim; ++i) - for (unsigned int j = 0; j < dim; ++j) - for (unsigned int l = 0; l < dim; ++l) - result[i][j][l] = (second[0][j][l] * - data.mapping_support_points[0][i]); - for (unsigned int k = 1; k < data.n_shape_functions; ++k) - for (unsigned int i = 0; i < spacedim; ++i) - for (unsigned int j = 0; j < dim; ++j) - for (unsigned int l = 0; l < dim; ++l) - result[i][j][l] += - (second[k][j][l] * - data.mapping_support_points[k][i]); - - // first push forward the j-components - for (unsigned int i = 0; i < spacedim; ++i) - for (unsigned int j = 0; j < spacedim; ++j) - for (unsigned int l = 0; l < dim; ++l) - { - tmp[i][j][l] = - result[i][0][l] * data.covariant[point][j][0]; - for (unsigned int jr = 1; jr < dim; ++jr) - { - tmp[i][j][l] += result[i][jr][l] * - data.covariant[point][j][jr]; - } - } - - // now, pushing forward the l-components - for (unsigned int i = 0; i < spacedim; ++i) - for (unsigned int j = 0; j < spacedim; ++j) - for (unsigned int l = 0; l < spacedim; ++l) - { - jacobian_pushed_forward_grads[point][i][j][l] = - tmp[i][j][0] * data.covariant[point][l][0]; - for (unsigned int lr = 1; lr < dim; ++lr) - { - jacobian_pushed_forward_grads[point][i][j][l] += - tmp[i][j][lr] * data.covariant[point][l][lr]; - } - } - } - } - } - } - - /** - * Update the third derivatives of the transformation from unit to real - * cell, the Jacobian hessians. - * - * Skip the computation if possible as indicated by the first argument. - */ - template - void - maybe_update_jacobian_2nd_derivatives( - const CellSimilarity::Similarity cell_similarity, - const typename QProjector::DataSetDescriptor data_set, - const typename dealii::MappingQGeneric::InternalData - & data, - std::vector> &jacobian_2nd_derivatives) - { - const UpdateFlags update_flags = data.update_each; - if (update_flags & update_jacobian_2nd_derivatives) - { - const unsigned int n_q_points = jacobian_2nd_derivatives.size(); - - if (cell_similarity != CellSimilarity::translation) - { - for (unsigned int point = 0; point < n_q_points; ++point) - { - const Tensor<3, dim> *third = - &data.third_derivative(point + data_set, 0); - double result[spacedim][dim][dim][dim]; - for (unsigned int i = 0; i < spacedim; ++i) - for (unsigned int j = 0; j < dim; ++j) - for (unsigned int l = 0; l < dim; ++l) - for (unsigned int m = 0; m < dim; ++m) - result[i][j][l][m] = - (third[0][j][l][m] * - data.mapping_support_points[0][i]); - for (unsigned int k = 1; k < data.n_shape_functions; ++k) - for (unsigned int i = 0; i < spacedim; ++i) - for (unsigned int j = 0; j < dim; ++j) - for (unsigned int l = 0; l < dim; ++l) - for (unsigned int m = 0; m < dim; ++m) - result[i][j][l][m] += - (third[k][j][l][m] * - data.mapping_support_points[k][i]); - - for (unsigned int i = 0; i < spacedim; ++i) - for (unsigned int j = 0; j < dim; ++j) - for (unsigned int l = 0; l < dim; ++l) - for (unsigned int m = 0; m < dim; ++m) - jacobian_2nd_derivatives[point][i][j][l][m] = - result[i][j][l][m]; - } - } - } - } - - /** - * Update the Hessian of the Hessian of the transformation from unit - * to real cell, the Jacobian Hessian gradients, pushed forward to the - * real cell coordinates. - * - * Skip the computation if possible as indicated by the first argument. - */ - template - void - maybe_update_jacobian_pushed_forward_2nd_derivatives( - const CellSimilarity::Similarity cell_similarity, - const typename QProjector::DataSetDescriptor data_set, - const typename dealii::MappingQGeneric::InternalData - &data, - std::vector> - &jacobian_pushed_forward_2nd_derivatives) - { - const UpdateFlags update_flags = data.update_each; - if (update_flags & update_jacobian_pushed_forward_2nd_derivatives) - { - const unsigned int n_q_points = - jacobian_pushed_forward_2nd_derivatives.size(); - - if (cell_similarity != CellSimilarity::translation) - { - double tmp[spacedim][spacedim][spacedim][spacedim]; - for (unsigned int point = 0; point < n_q_points; ++point) - { - const Tensor<3, dim> *third = - &data.third_derivative(point + data_set, 0); - double result[spacedim][dim][dim][dim]; - for (unsigned int i = 0; i < spacedim; ++i) - for (unsigned int j = 0; j < dim; ++j) - for (unsigned int l = 0; l < dim; ++l) - for (unsigned int m = 0; m < dim; ++m) - result[i][j][l][m] = - (third[0][j][l][m] * - data.mapping_support_points[0][i]); - for (unsigned int k = 1; k < data.n_shape_functions; ++k) - for (unsigned int i = 0; i < spacedim; ++i) - for (unsigned int j = 0; j < dim; ++j) - for (unsigned int l = 0; l < dim; ++l) - for (unsigned int m = 0; m < dim; ++m) - result[i][j][l][m] += - (third[k][j][l][m] * - data.mapping_support_points[k][i]); - - // push forward the j-coordinate - for (unsigned int i = 0; i < spacedim; ++i) - for (unsigned int j = 0; j < spacedim; ++j) - for (unsigned int l = 0; l < dim; ++l) - for (unsigned int m = 0; m < dim; ++m) - { - jacobian_pushed_forward_2nd_derivatives - [point][i][j][l][m] = - result[i][0][l][m] * - data.covariant[point][j][0]; - for (unsigned int jr = 1; jr < dim; ++jr) - jacobian_pushed_forward_2nd_derivatives[point] - [i][j][l] - [m] += - result[i][jr][l][m] * - data.covariant[point][j][jr]; - } - - // push forward the l-coordinate - for (unsigned int i = 0; i < spacedim; ++i) - for (unsigned int j = 0; j < spacedim; ++j) - for (unsigned int l = 0; l < spacedim; ++l) - for (unsigned int m = 0; m < dim; ++m) - { - tmp[i][j][l][m] = - jacobian_pushed_forward_2nd_derivatives[point] - [i][j][0] - [m] * - data.covariant[point][l][0]; - for (unsigned int lr = 1; lr < dim; ++lr) - tmp[i][j][l][m] += - jacobian_pushed_forward_2nd_derivatives - [point][i][j][lr][m] * - data.covariant[point][l][lr]; - } - - // push forward the m-coordinate - for (unsigned int i = 0; i < spacedim; ++i) - for (unsigned int j = 0; j < spacedim; ++j) - for (unsigned int l = 0; l < spacedim; ++l) - for (unsigned int m = 0; m < spacedim; ++m) - { - jacobian_pushed_forward_2nd_derivatives - [point][i][j][l][m] = - tmp[i][j][l][0] * data.covariant[point][m][0]; - for (unsigned int mr = 1; mr < dim; ++mr) - jacobian_pushed_forward_2nd_derivatives[point] - [i][j][l] - [m] += - tmp[i][j][l][mr] * - data.covariant[point][m][mr]; - } - } - } - } - } - - /** - * Update the fourth derivatives of the transformation from unit to real - * cell, the Jacobian hessian gradients. - * - * Skip the computation if possible as indicated by the first argument. - */ - template - void - maybe_update_jacobian_3rd_derivatives( - const CellSimilarity::Similarity cell_similarity, - const typename QProjector::DataSetDescriptor data_set, - const typename dealii::MappingQGeneric::InternalData - & data, - std::vector> &jacobian_3rd_derivatives) - { - const UpdateFlags update_flags = data.update_each; - if (update_flags & update_jacobian_3rd_derivatives) - { - const unsigned int n_q_points = jacobian_3rd_derivatives.size(); - - if (cell_similarity != CellSimilarity::translation) - { - for (unsigned int point = 0; point < n_q_points; ++point) - { - const Tensor<4, dim> *fourth = - &data.fourth_derivative(point + data_set, 0); - double result[spacedim][dim][dim][dim][dim]; - for (unsigned int i = 0; i < spacedim; ++i) - for (unsigned int j = 0; j < dim; ++j) - for (unsigned int l = 0; l < dim; ++l) - for (unsigned int m = 0; m < dim; ++m) - for (unsigned int n = 0; n < dim; ++n) - result[i][j][l][m][n] = - (fourth[0][j][l][m][n] * - data.mapping_support_points[0][i]); - for (unsigned int k = 1; k < data.n_shape_functions; ++k) - for (unsigned int i = 0; i < spacedim; ++i) - for (unsigned int j = 0; j < dim; ++j) - for (unsigned int l = 0; l < dim; ++l) - for (unsigned int m = 0; m < dim; ++m) - for (unsigned int n = 0; n < dim; ++n) - result[i][j][l][m][n] += - (fourth[k][j][l][m][n] * - data.mapping_support_points[k][i]); - - for (unsigned int i = 0; i < spacedim; ++i) - for (unsigned int j = 0; j < dim; ++j) - for (unsigned int l = 0; l < dim; ++l) - for (unsigned int m = 0; m < dim; ++m) - for (unsigned int n = 0; n < dim; ++n) - jacobian_3rd_derivatives[point][i][j][l][m][n] = - result[i][j][l][m][n]; - } - } - } - } - - /** - * Update the Hessian gradient of the transformation from unit to real - * cell, the Jacobian Hessians, pushed forward to the real cell - * coordinates. - * - * Skip the computation if possible as indicated by the first argument. - */ - template - void - maybe_update_jacobian_pushed_forward_3rd_derivatives( - const CellSimilarity::Similarity cell_similarity, - const typename QProjector::DataSetDescriptor data_set, - const typename dealii::MappingQGeneric::InternalData - &data, - std::vector> - &jacobian_pushed_forward_3rd_derivatives) - { - const UpdateFlags update_flags = data.update_each; - if (update_flags & update_jacobian_pushed_forward_3rd_derivatives) - { - const unsigned int n_q_points = - jacobian_pushed_forward_3rd_derivatives.size(); - - if (cell_similarity != CellSimilarity::translation) - { - double tmp[spacedim][spacedim][spacedim][spacedim][spacedim]; - for (unsigned int point = 0; point < n_q_points; ++point) - { - const Tensor<4, dim> *fourth = - &data.fourth_derivative(point + data_set, 0); - double result[spacedim][dim][dim][dim][dim]; - for (unsigned int i = 0; i < spacedim; ++i) - for (unsigned int j = 0; j < dim; ++j) - for (unsigned int l = 0; l < dim; ++l) - for (unsigned int m = 0; m < dim; ++m) - for (unsigned int n = 0; n < dim; ++n) - result[i][j][l][m][n] = - (fourth[0][j][l][m][n] * - data.mapping_support_points[0][i]); - for (unsigned int k = 1; k < data.n_shape_functions; ++k) - for (unsigned int i = 0; i < spacedim; ++i) - for (unsigned int j = 0; j < dim; ++j) - for (unsigned int l = 0; l < dim; ++l) - for (unsigned int m = 0; m < dim; ++m) - for (unsigned int n = 0; n < dim; ++n) - result[i][j][l][m][n] += - (fourth[k][j][l][m][n] * - data.mapping_support_points[k][i]); - - // push-forward the j-coordinate - for (unsigned int i = 0; i < spacedim; ++i) - for (unsigned int j = 0; j < spacedim; ++j) - for (unsigned int l = 0; l < dim; ++l) - for (unsigned int m = 0; m < dim; ++m) - for (unsigned int n = 0; n < dim; ++n) - { - tmp[i][j][l][m][n] = - result[i][0][l][m][n] * - data.covariant[point][j][0]; - for (unsigned int jr = 1; jr < dim; ++jr) - tmp[i][j][l][m][n] += - result[i][jr][l][m][n] * - data.covariant[point][j][jr]; - } - - // push-forward the l-coordinate - for (unsigned int i = 0; i < spacedim; ++i) - for (unsigned int j = 0; j < spacedim; ++j) - for (unsigned int l = 0; l < spacedim; ++l) - for (unsigned int m = 0; m < dim; ++m) - for (unsigned int n = 0; n < dim; ++n) - { - jacobian_pushed_forward_3rd_derivatives - [point][i][j][l][m][n] = - tmp[i][j][0][m][n] * - data.covariant[point][l][0]; - for (unsigned int lr = 1; lr < dim; ++lr) - jacobian_pushed_forward_3rd_derivatives - [point][i][j][l][m][n] += - tmp[i][j][lr][m][n] * - data.covariant[point][l][lr]; - } - - // push-forward the m-coordinate - for (unsigned int i = 0; i < spacedim; ++i) - for (unsigned int j = 0; j < spacedim; ++j) - for (unsigned int l = 0; l < spacedim; ++l) - for (unsigned int m = 0; m < spacedim; ++m) - for (unsigned int n = 0; n < dim; ++n) - { - tmp[i][j][l][m][n] = - jacobian_pushed_forward_3rd_derivatives - [point][i][j][l][0][n] * - data.covariant[point][m][0]; - for (unsigned int mr = 1; mr < dim; ++mr) - tmp[i][j][l][m][n] += - jacobian_pushed_forward_3rd_derivatives - [point][i][j][l][mr][n] * - data.covariant[point][m][mr]; - } - - // push-forward the n-coordinate - for (unsigned int i = 0; i < spacedim; ++i) - for (unsigned int j = 0; j < spacedim; ++j) - for (unsigned int l = 0; l < spacedim; ++l) - for (unsigned int m = 0; m < spacedim; ++m) - for (unsigned int n = 0; n < spacedim; ++n) - { - jacobian_pushed_forward_3rd_derivatives - [point][i][j][l][m][n] = - tmp[i][j][l][m][0] * - data.covariant[point][n][0]; - for (unsigned int nr = 1; nr < dim; ++nr) - jacobian_pushed_forward_3rd_derivatives - [point][i][j][l][m][n] += - tmp[i][j][l][m][nr] * - data.covariant[point][n][nr]; - } - } - } - } - } - } // namespace - } // namespace MappingQGenericImplementation -} // namespace internal - - - template MappingQGeneric::MappingQGeneric(const unsigned int p) : polynomial_degree(p) @@ -2674,265 +1135,6 @@ MappingQGeneric::fill_fe_values( -namespace internal -{ - namespace MappingQGenericImplementation - { - namespace - { - /** - * Depending on what information is called for in the update flags of the - * @p data object, compute the various pieces of information that is - * required by the fill_fe_face_values() and fill_fe_subface_values() - * functions. This function simply unifies the work that would be done by - * those two functions. - * - * The resulting data is put into the @p output_data argument. - */ - template - void - maybe_compute_face_data( - const dealii::MappingQGeneric &mapping, - const typename dealii::Triangulation::cell_iterator - & cell, - const unsigned int face_no, - const unsigned int subface_no, - const unsigned int n_q_points, - const std::vector &weights, - const typename dealii::MappingQGeneric::InternalData - &data, - internal::FEValuesImplementation::MappingRelatedData - &output_data) - { - const UpdateFlags update_flags = data.update_each; - - if (update_flags & - (update_boundary_forms | update_normal_vectors | update_jacobians | - update_JxW_values | update_inverse_jacobians)) - { - if (update_flags & update_boundary_forms) - AssertDimension(output_data.boundary_forms.size(), n_q_points); - if (update_flags & update_normal_vectors) - AssertDimension(output_data.normal_vectors.size(), n_q_points); - if (update_flags & update_JxW_values) - AssertDimension(output_data.JxW_values.size(), n_q_points); - - Assert(data.aux.size() + 1 >= dim, ExcInternalError()); - - // first compute some common data that is used for evaluating - // all of the flags below - - // map the unit tangentials to the real cell. checking for d!=dim-1 - // eliminates compiler warnings regarding unsigned int expressions < - // 0. - for (unsigned int d = 0; d != dim - 1; ++d) - { - Assert(face_no + GeometryInfo::faces_per_cell * d < - data.unit_tangentials.size(), - ExcInternalError()); - Assert( - data.aux[d].size() <= - data - .unit_tangentials[face_no + - GeometryInfo::faces_per_cell * d] - .size(), - ExcInternalError()); - - mapping.transform( - make_array_view( - data - .unit_tangentials[face_no + - GeometryInfo::faces_per_cell * d]), - mapping_contravariant, - data, - make_array_view(data.aux[d])); - } - - if (update_flags & update_boundary_forms) - { - // if dim==spacedim, we can use the unit tangentials to compute - // the boundary form by simply taking the cross product - if (dim == spacedim) - { - for (unsigned int i = 0; i < n_q_points; ++i) - switch (dim) - { - case 1: - // in 1d, we don't have access to any of the - // data.aux fields (because it has only dim-1 - // components), but we can still compute the - // boundary form by simply looking at the number of - // the face - output_data.boundary_forms[i][0] = - (face_no == 0 ? -1 : +1); - break; - case 2: - output_data.boundary_forms[i] = - cross_product_2d(data.aux[0][i]); - break; - case 3: - output_data.boundary_forms[i] = - cross_product_3d(data.aux[0][i], data.aux[1][i]); - break; - default: - Assert(false, ExcNotImplemented()); - } - } - else //(dim < spacedim) - { - // in the codim-one case, the boundary form results from the - // cross product of all the face tangential vectors and the - // cell normal vector - // - // to compute the cell normal, use the same method used in - // fill_fe_values for cells above - AssertDimension(data.contravariant.size(), n_q_points); - - for (unsigned int point = 0; point < n_q_points; ++point) - { - if (dim == 1) - { - // J is a tangent vector - output_data.boundary_forms[point] = - data.contravariant[point].transpose()[0]; - output_data.boundary_forms[point] /= - (face_no == 0 ? -1. : +1.) * - output_data.boundary_forms[point].norm(); - } - - if (dim == 2) - { - const DerivativeForm<1, spacedim, dim> DX_t = - data.contravariant[point].transpose(); - - Tensor<1, spacedim> cell_normal = - cross_product_3d(DX_t[0], DX_t[1]); - cell_normal /= cell_normal.norm(); - - // then compute the face normal from the face - // tangent and the cell normal: - output_data.boundary_forms[point] = - cross_product_3d(data.aux[0][point], cell_normal); - } - } - } - } - - if (update_flags & update_JxW_values) - for (unsigned int i = 0; i < output_data.boundary_forms.size(); - ++i) - { - output_data.JxW_values[i] = - output_data.boundary_forms[i].norm() * weights[i]; - - if (subface_no != numbers::invalid_unsigned_int) - { - const double area_ratio = - GeometryInfo::subface_ratio( - cell->subface_case(face_no), subface_no); - output_data.JxW_values[i] *= area_ratio; - } - } - - if (update_flags & update_normal_vectors) - for (unsigned int i = 0; i < output_data.normal_vectors.size(); - ++i) - output_data.normal_vectors[i] = - Point(output_data.boundary_forms[i] / - output_data.boundary_forms[i].norm()); - - if (update_flags & update_jacobians) - for (unsigned int point = 0; point < n_q_points; ++point) - output_data.jacobians[point] = data.contravariant[point]; - - if (update_flags & update_inverse_jacobians) - for (unsigned int point = 0; point < n_q_points; ++point) - output_data.inverse_jacobians[point] = - data.covariant[point].transpose(); - } - } - - - /** - * Do the work of MappingQGeneric::fill_fe_face_values() and - * MappingQGeneric::fill_fe_subface_values() in a generic way, - * using the 'data_set' to differentiate whether we will - * work on a face (and if so, which one) or subface. - */ - template - void - do_fill_fe_face_values( - const dealii::MappingQGeneric &mapping, - const typename dealii::Triangulation::cell_iterator - & cell, - const unsigned int face_no, - const unsigned int subface_no, - const typename QProjector::DataSetDescriptor data_set, - const Quadrature & quadrature, - const typename dealii::MappingQGeneric::InternalData - &data, - internal::FEValuesImplementation::MappingRelatedData - &output_data) - { - if (dim > 1 && data.tensor_product_quadrature) - { - maybe_update_q_points_Jacobians_and_grads_tensor( - CellSimilarity::none, - data, - output_data.quadrature_points, - output_data.jacobian_grads); - } - else - { - maybe_compute_q_points( - data_set, data, output_data.quadrature_points); - maybe_update_Jacobians(CellSimilarity::none, - data_set, - data); - maybe_update_jacobian_grads( - CellSimilarity::none, data_set, data, output_data.jacobian_grads); - } - maybe_update_jacobian_pushed_forward_grads( - CellSimilarity::none, - data_set, - data, - output_data.jacobian_pushed_forward_grads); - maybe_update_jacobian_2nd_derivatives( - CellSimilarity::none, - data_set, - data, - output_data.jacobian_2nd_derivatives); - maybe_update_jacobian_pushed_forward_2nd_derivatives( - CellSimilarity::none, - data_set, - data, - output_data.jacobian_pushed_forward_2nd_derivatives); - maybe_update_jacobian_3rd_derivatives( - CellSimilarity::none, - data_set, - data, - output_data.jacobian_3rd_derivatives); - maybe_update_jacobian_pushed_forward_3rd_derivatives( - CellSimilarity::none, - data_set, - data, - output_data.jacobian_pushed_forward_3rd_derivatives); - - maybe_compute_face_data(mapping, - cell, - face_no, - subface_no, - quadrature.size(), - quadrature.get_weights(), - data, - output_data); - } - } // namespace - } // namespace MappingQGenericImplementation -} // namespace internal - - - template void MappingQGeneric::fill_fe_face_values( @@ -3029,403 +1231,6 @@ MappingQGeneric::fill_fe_subface_values( -namespace internal -{ - namespace MappingQGenericImplementation - { - namespace - { - template - void - transform_fields( - const ArrayView> & input, - const MappingKind mapping_kind, - const typename Mapping::InternalDataBase &mapping_data, - const ArrayView> & output) - { - AssertDimension(input.size(), output.size()); - Assert((dynamic_cast::InternalData *>( - &mapping_data) != nullptr), - ExcInternalError()); - const typename dealii::MappingQGeneric::InternalData - &data = - static_cast:: - InternalData &>(mapping_data); - - switch (mapping_kind) - { - case mapping_contravariant: - { - Assert( - data.update_each & update_contravariant_transformation, - typename FEValuesBase::ExcAccessToUninitializedField( - "update_contravariant_transformation")); - - for (unsigned int i = 0; i < output.size(); ++i) - output[i] = - apply_transformation(data.contravariant[i], input[i]); - - return; - } - - case mapping_piola: - { - Assert( - data.update_each & update_contravariant_transformation, - typename FEValuesBase::ExcAccessToUninitializedField( - "update_contravariant_transformation")); - Assert( - data.update_each & update_volume_elements, - typename FEValuesBase::ExcAccessToUninitializedField( - "update_volume_elements")); - Assert(rank == 1, ExcMessage("Only for rank 1")); - if (rank != 1) - return; - - for (unsigned int i = 0; i < output.size(); ++i) - { - output[i] = - apply_transformation(data.contravariant[i], input[i]); - output[i] /= data.volume_elements[i]; - } - return; - } - // We still allow this operation as in the - // reference cell Derivatives are Tensor - // rather than DerivativeForm - case mapping_covariant: - { - Assert( - data.update_each & update_contravariant_transformation, - typename FEValuesBase::ExcAccessToUninitializedField( - "update_covariant_transformation")); - - for (unsigned int i = 0; i < output.size(); ++i) - output[i] = apply_transformation(data.covariant[i], input[i]); - - return; - } - - default: - Assert(false, ExcNotImplemented()); - } - } - - - template - void - transform_gradients( - const ArrayView> & input, - const MappingKind mapping_kind, - const typename Mapping::InternalDataBase &mapping_data, - const ArrayView> & output) - { - AssertDimension(input.size(), output.size()); - Assert((dynamic_cast::InternalData *>( - &mapping_data) != nullptr), - ExcInternalError()); - const typename dealii::MappingQGeneric::InternalData - &data = - static_cast:: - InternalData &>(mapping_data); - - switch (mapping_kind) - { - case mapping_contravariant_gradient: - { - Assert( - data.update_each & update_covariant_transformation, - typename FEValuesBase::ExcAccessToUninitializedField( - "update_covariant_transformation")); - Assert( - data.update_each & update_contravariant_transformation, - typename FEValuesBase::ExcAccessToUninitializedField( - "update_contravariant_transformation")); - Assert(rank == 2, ExcMessage("Only for rank 2")); - - for (unsigned int i = 0; i < output.size(); ++i) - { - const DerivativeForm<1, spacedim, dim> A = - apply_transformation(data.contravariant[i], - transpose(input[i])); - output[i] = - apply_transformation(data.covariant[i], A.transpose()); - } - - return; - } - - case mapping_covariant_gradient: - { - Assert( - data.update_each & update_covariant_transformation, - typename FEValuesBase::ExcAccessToUninitializedField( - "update_covariant_transformation")); - Assert(rank == 2, ExcMessage("Only for rank 2")); - - for (unsigned int i = 0; i < output.size(); ++i) - { - const DerivativeForm<1, spacedim, dim> A = - apply_transformation(data.covariant[i], - transpose(input[i])); - output[i] = - apply_transformation(data.covariant[i], A.transpose()); - } - - return; - } - - case mapping_piola_gradient: - { - Assert( - data.update_each & update_covariant_transformation, - typename FEValuesBase::ExcAccessToUninitializedField( - "update_covariant_transformation")); - Assert( - data.update_each & update_contravariant_transformation, - typename FEValuesBase::ExcAccessToUninitializedField( - "update_contravariant_transformation")); - Assert( - data.update_each & update_volume_elements, - typename FEValuesBase::ExcAccessToUninitializedField( - "update_volume_elements")); - Assert(rank == 2, ExcMessage("Only for rank 2")); - - for (unsigned int i = 0; i < output.size(); ++i) - { - const DerivativeForm<1, spacedim, dim> A = - apply_transformation(data.covariant[i], input[i]); - const Tensor<2, spacedim> T = - apply_transformation(data.contravariant[i], - A.transpose()); - - output[i] = transpose(T); - output[i] /= data.volume_elements[i]; - } - - return; - } - - default: - Assert(false, ExcNotImplemented()); - } - } - - - - template - void - transform_hessians( - const ArrayView> & input, - const MappingKind mapping_kind, - const typename Mapping::InternalDataBase &mapping_data, - const ArrayView> & output) - { - AssertDimension(input.size(), output.size()); - Assert((dynamic_cast::InternalData *>( - &mapping_data) != nullptr), - ExcInternalError()); - const typename dealii::MappingQGeneric::InternalData - &data = - static_cast:: - InternalData &>(mapping_data); - - switch (mapping_kind) - { - case mapping_contravariant_hessian: - { - Assert( - data.update_each & update_covariant_transformation, - typename FEValuesBase::ExcAccessToUninitializedField( - "update_covariant_transformation")); - Assert( - data.update_each & update_contravariant_transformation, - typename FEValuesBase::ExcAccessToUninitializedField( - "update_contravariant_transformation")); - - for (unsigned int q = 0; q < output.size(); ++q) - for (unsigned int i = 0; i < spacedim; ++i) - { - double tmp1[dim][dim]; - for (unsigned int J = 0; J < dim; ++J) - for (unsigned int K = 0; K < dim; ++K) - { - tmp1[J][K] = - data.contravariant[q][i][0] * input[q][0][J][K]; - for (unsigned int I = 1; I < dim; ++I) - tmp1[J][K] += - data.contravariant[q][i][I] * input[q][I][J][K]; - } - for (unsigned int j = 0; j < spacedim; ++j) - { - double tmp2[dim]; - for (unsigned int K = 0; K < dim; ++K) - { - tmp2[K] = data.covariant[q][j][0] * tmp1[0][K]; - for (unsigned int J = 1; J < dim; ++J) - tmp2[K] += data.covariant[q][j][J] * tmp1[J][K]; - } - for (unsigned int k = 0; k < spacedim; ++k) - { - output[q][i][j][k] = - data.covariant[q][k][0] * tmp2[0]; - for (unsigned int K = 1; K < dim; ++K) - output[q][i][j][k] += - data.covariant[q][k][K] * tmp2[K]; - } - } - } - return; - } - - case mapping_covariant_hessian: - { - Assert( - data.update_each & update_covariant_transformation, - typename FEValuesBase::ExcAccessToUninitializedField( - "update_covariant_transformation")); - - for (unsigned int q = 0; q < output.size(); ++q) - for (unsigned int i = 0; i < spacedim; ++i) - { - double tmp1[dim][dim]; - for (unsigned int J = 0; J < dim; ++J) - for (unsigned int K = 0; K < dim; ++K) - { - tmp1[J][K] = - data.covariant[q][i][0] * input[q][0][J][K]; - for (unsigned int I = 1; I < dim; ++I) - tmp1[J][K] += - data.covariant[q][i][I] * input[q][I][J][K]; - } - for (unsigned int j = 0; j < spacedim; ++j) - { - double tmp2[dim]; - for (unsigned int K = 0; K < dim; ++K) - { - tmp2[K] = data.covariant[q][j][0] * tmp1[0][K]; - for (unsigned int J = 1; J < dim; ++J) - tmp2[K] += data.covariant[q][j][J] * tmp1[J][K]; - } - for (unsigned int k = 0; k < spacedim; ++k) - { - output[q][i][j][k] = - data.covariant[q][k][0] * tmp2[0]; - for (unsigned int K = 1; K < dim; ++K) - output[q][i][j][k] += - data.covariant[q][k][K] * tmp2[K]; - } - } - } - - return; - } - - case mapping_piola_hessian: - { - Assert( - data.update_each & update_covariant_transformation, - typename FEValuesBase::ExcAccessToUninitializedField( - "update_covariant_transformation")); - Assert( - data.update_each & update_contravariant_transformation, - typename FEValuesBase::ExcAccessToUninitializedField( - "update_contravariant_transformation")); - Assert( - data.update_each & update_volume_elements, - typename FEValuesBase::ExcAccessToUninitializedField( - "update_volume_elements")); - - for (unsigned int q = 0; q < output.size(); ++q) - for (unsigned int i = 0; i < spacedim; ++i) - { - double factor[dim]; - for (unsigned int I = 0; I < dim; ++I) - factor[I] = - data.contravariant[q][i][I] / data.volume_elements[q]; - double tmp1[dim][dim]; - for (unsigned int J = 0; J < dim; ++J) - for (unsigned int K = 0; K < dim; ++K) - { - tmp1[J][K] = factor[0] * input[q][0][J][K]; - for (unsigned int I = 1; I < dim; ++I) - tmp1[J][K] += factor[I] * input[q][I][J][K]; - } - for (unsigned int j = 0; j < spacedim; ++j) - { - double tmp2[dim]; - for (unsigned int K = 0; K < dim; ++K) - { - tmp2[K] = data.covariant[q][j][0] * tmp1[0][K]; - for (unsigned int J = 1; J < dim; ++J) - tmp2[K] += data.covariant[q][j][J] * tmp1[J][K]; - } - for (unsigned int k = 0; k < spacedim; ++k) - { - output[q][i][j][k] = - data.covariant[q][k][0] * tmp2[0]; - for (unsigned int K = 1; K < dim; ++K) - output[q][i][j][k] += - data.covariant[q][k][K] * tmp2[K]; - } - } - } - - return; - } - - default: - Assert(false, ExcNotImplemented()); - } - } - - - - template - void - transform_differential_forms( - const ArrayView> &input, - const MappingKind mapping_kind, - const typename Mapping::InternalDataBase &mapping_data, - const ArrayView> & output) - { - AssertDimension(input.size(), output.size()); - Assert((dynamic_cast::InternalData *>( - &mapping_data) != nullptr), - ExcInternalError()); - const typename dealii::MappingQGeneric::InternalData - &data = - static_cast:: - InternalData &>(mapping_data); - - switch (mapping_kind) - { - case mapping_covariant: - { - Assert( - data.update_each & update_contravariant_transformation, - typename FEValuesBase::ExcAccessToUninitializedField( - "update_covariant_transformation")); - - for (unsigned int i = 0; i < output.size(); ++i) - output[i] = apply_transformation(data.covariant[i], input[i]); - - return; - } - default: - Assert(false, ExcNotImplemented()); - } - } - } // namespace - } // namespace MappingQGenericImplementation -} // namespace internal - - - template void MappingQGeneric::transform( -- 2.39.5