From: Martin Kronbichler Date: Fri, 13 Mar 2020 08:27:21 +0000 (+0100) Subject: Implement matrix-free evaluation of geometry for MF::MappingInfo X-Git-Tag: v9.2.0-rc1~342^2~2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=4606f073f13d3f138b3904c2addfb24326200945;p=dealii.git Implement matrix-free evaluation of geometry for MF::MappingInfo --- diff --git a/include/deal.II/matrix_free/fe_evaluation.h b/include/deal.II/matrix_free/fe_evaluation.h index c9d6778df3..71bd3457df 100644 --- a/include/deal.II/matrix_free/fe_evaluation.h +++ b/include/deal.II/matrix_free/fe_evaluation.h @@ -6928,39 +6928,34 @@ FEEvaluationdata->data.front().n_q_points_1d : n_q_points_1d; - - // Cartesian mesh: not all quadrature points are stored, only the - // diagonal. Hence, need to find the tensor product index and retrieve the - // value from that const Point *quadrature_points = &this->mapping_data->quadrature_points [this->mapping_data->quadrature_point_offsets[this->cell]]; - if (this->cell_type == internal::MatrixFreeFunctions::cartesian) + + // Cartesian/affine mesh: only first vertex of cell is stored, we must + // compute it through the Jacobian (which is stored in non-inverted and + // non-transposed form as index '1' in the jacobian field) + if (this->cell_type <= internal::MatrixFreeFunctions::affine) { - Point point; - switch (dim) - { - case 1: - return quadrature_points[q]; - case 2: - point[0] = quadrature_points[q % n_q_points_1d_actual][0]; - point[1] = quadrature_points[q / n_q_points_1d_actual][1]; - return point; - case 3: - point[0] = quadrature_points[q % n_q_points_1d_actual][0]; - point[1] = quadrature_points[(q / n_q_points_1d_actual) % - n_q_points_1d_actual][1]; - point[2] = quadrature_points[q / (n_q_points_1d_actual * - n_q_points_1d_actual)][2]; - return point; - default: - Assert(false, ExcNotImplemented()); - return point; - } + Assert(this->jacobian != nullptr, ExcNotInitialized()); + Point point = quadrature_points[0]; + + const Tensor<2, dim, VectorizedArrayType> &jac = this->jacobian[1]; + if (this->cell_type == internal::MatrixFreeFunctions::cartesian) + for (unsigned int d = 0; d < dim; ++d) + point[d] += jac[d][d] * + static_cast( + this->mapping_data->descriptor[this->active_quad_index] + .quadrature.point(q)[d]); + else + for (unsigned int d = 0; d < dim; ++d) + for (unsigned int e = 0; e < dim; ++e) + point[d] += jac[d][e] * static_cast( + this->mapping_data + ->descriptor[this->active_quad_index] + .quadrature.point(q)[e]); + return point; } - // all other cases: just return the respective data as it is fully stored else return quadrature_points[q]; } diff --git a/include/deal.II/matrix_free/mapping_info.h b/include/deal.II/matrix_free/mapping_info.h index 120fa1379c..30040fd4a7 100644 --- a/include/deal.II/matrix_free/mapping_info.h +++ b/include/deal.II/matrix_free/mapping_info.h @@ -441,6 +441,34 @@ namespace internal */ SmartPointer> mapping; + /** + * Internal function to compute the geometry for the case the mapping is + * a MappingQ and a single quadrature formula per slot (non-hp case) is + * used. This method computes all data from the underlying cell + * quadrature points using the fast operator evaluation techniques from + * the matrix-free framework itself, i.e., it uses a polynomial + * description of the cell geometry (that is computed in a first step) + * and then computes all Jacobians and normal vectors based on this + * information. This optimized approach is much faster than going + * through FEValues and FEFaceValues, especially when several different + * quadrature formulas are involved, and consumes less memory. + * + * @param tria The triangulation to be used for setup + * + * @param cells The actual cells of the triangulation to be worked on, + * given as a tuple of the level and index within the level as used in + * the main initialization of the class + * + * @param faces The description of the connectivity from faces to cells + * as filled in the MatrixFree class + */ + void + compute_mapping_q( + const dealii::Triangulation & tria, + const std::vector> &cells, + const std::vector< + FaceToCellTopology> &faces); + /** * Computes the information in the given cells, called within * initialize. @@ -529,7 +557,7 @@ namespace internal * comparator class within a std::map<> of the given arrays. Note that this * comparison operator does not satisfy all the mathematical properties one * usually wants to have (consider e.g. the numbers a=0, b=0.1, c=0.2 with - * tolerance 0.15; the operator gives a &v1, const std::vector &v2) const; + /** + * Compare two vectorized arrays (stored as tensors to avoid alignment + * issues). + */ bool operator()( const Tensor<1, VectorizedArrayType::size(), Number> &t1, const Tensor<1, VectorizedArrayType::size(), Number> &t2) const; + /** + * Compare two rank-1 tensors of vectorized arrays (stored as tensors to + * avoid alignment issues). + */ template bool operator()( @@ -556,6 +595,10 @@ namespace internal const Tensor<1, dim, Tensor<1, VectorizedArrayType::size(), Number>> &t2) const; + /** + * Compare two rank-2 tensors of vectorized arrays (stored as tensors to + * avoid alignment issues). + */ template bool operator()( @@ -564,6 +607,14 @@ namespace internal const Tensor<2, dim, Tensor<1, VectorizedArrayType::size(), Number>> &t2) const; + /** + * Compare two arrays of tensors. + */ + template + bool + operator()(const std::array, dim + 1> &t1, + const std::array, dim + 1> &t2) const; + Number tolerance; }; diff --git a/include/deal.II/matrix_free/mapping_info.templates.h b/include/deal.II/matrix_free/mapping_info.templates.h index bb41fc52ac..5788897420 100644 --- a/include/deal.II/matrix_free/mapping_info.templates.h +++ b/include/deal.II/matrix_free/mapping_info.templates.h @@ -23,13 +23,15 @@ #include #include +#include #include #include -#include +#include +#include +#include #include - DEAL_II_NAMESPACE_OPEN @@ -243,6 +245,7 @@ namespace internal face_data_by_cells.clear(); cell_type.clear(); face_type.clear(); + mapping = nullptr; } @@ -349,11 +352,20 @@ namespace internal update_default); } - // Could call these functions in parallel, but not useful because the - // work inside is nicely split up already - initialize_cells(tria, cells, active_fe_index, mapping); - initialize_faces(tria, cells, face_info.faces, mapping); - initialize_faces_by_cells(tria, cells, mapping); + // In case we have no hp adaptivity (active_fe_index is empty), we have + // cells, and the mapping is MappingQGeneric or a derived class, we can + // use the fast method. + if (active_fe_index.empty() && !cells.empty() && + dynamic_cast *>(&mapping)) + compute_mapping_q(tria, cells, face_info.faces); + else + { + // Could call these functions in parallel, but not useful because + // the work inside is nicely split up already + initialize_cells(tria, cells, active_fe_index, mapping); + initialize_faces(tria, cells, face_info.faces, mapping); + initialize_faces_by_cells(tria, cells, mapping); + } } @@ -379,17 +391,41 @@ namespace internal this->mapping = &mapping; - // Could call these functions in parallel, but not useful because the - // work inside is nicely split up already - initialize_cells(tria, cells, active_fe_index, mapping); - initialize_faces(tria, cells, face_info.faces, mapping); - initialize_faces_by_cells(tria, cells, mapping); + if (active_fe_index.empty() && !cells.empty() && + dynamic_cast *>(&mapping)) + compute_mapping_q(tria, cells, face_info.faces); + else + { + // Could call these functions in parallel, but not useful because + // the work inside is nicely split up already + initialize_cells(tria, cells, active_fe_index, mapping); + initialize_faces(tria, cells, face_info.faces, mapping); + initialize_faces_by_cells(tria, cells, mapping); + } } /* ------------------------- initialization of cells ------------------- */ + // Copy a vectorized array of one type to another type + template + inline DEAL_II_ALWAYS_INLINE void + store_vectorized_array(const VectorizedArrayType1 value, + const unsigned int offset, + VectorizedArrayType2 & result) + { + static_assert(VectorizedArrayType2::n_array_elements >= + VectorizedArrayType1::n_array_elements, + "Cannot convert to vectorized array of wider number type"); + + DEAL_II_OPENMP_SIMD_PRAGMA + for (unsigned int v = 0; v < VectorizedArrayType1::n_array_elements; ++v) + result[offset + v] = value[v]; + } + + + // Namespace with implementation of extraction of values on cell // range namespace ExtractCellHelper @@ -463,6 +499,69 @@ namespace internal } } + // For second derivatives on the real cell, we need the gradient of the + // inverse Jacobian J. This involves some calculus and is done + // vectorized. If L is the gradient of the jacobian on the unit cell, + // the gradient of the inverse is given by (multidimensional calculus) - + // J * (J * L) * J (the third J is because we need to transform the + // gradient L from the unit to the real cell, and then apply the inverse + // Jacobian). Compare this with 1D with j(x) = 1/k(phi(x)), where j = + // phi' is the inverse of the jacobian and k is the derivative of the + // jacobian on the unit cell. Then j' = phi' k'/k^2 = j k' j^2. + template + Tensor<1, dim *(dim + 1) / 2, Tensor<1, dim, Number>> + process_jacobian_gradient(const Tensor<2, dim, Number> &inv_jac, + const Tensor<3, dim, Number> &jac_grad) + { + Number inv_jac_grad[dim][dim][dim]; + + // compute: inv_jac_grad = J*grad_unit(J^-1) + for (unsigned int d = 0; d < dim; ++d) + for (unsigned int e = 0; e < dim; ++e) + for (unsigned int f = 0; f < dim; ++f) + { + inv_jac_grad[f][e][d] = (inv_jac[f][0] * jac_grad[d][e][0]); + for (unsigned int g = 1; g < dim; ++g) + inv_jac_grad[f][e][d] += (inv_jac[f][g] * jac_grad[d][e][g]); + } + + // compute: transpose (-jac * jac_grad[d] * jac) + Number tmp[dim]; + Number grad_jac_inv[dim][dim][dim]; + for (unsigned int d = 0; d < dim; ++d) + for (unsigned int e = 0; e < dim; ++e) + { + for (unsigned int f = 0; f < dim; ++f) + { + tmp[f] = Number(); + for (unsigned int g = 0; g < dim; ++g) + tmp[f] -= inv_jac_grad[d][f][g] * inv_jac[g][e]; + } + + // needed for non-diagonal part of Jacobian grad + for (unsigned int f = 0; f < dim; ++f) + { + grad_jac_inv[f][d][e] = inv_jac[f][0] * tmp[0]; + for (unsigned int g = 1; g < dim; ++g) + grad_jac_inv[f][d][e] += inv_jac[f][g] * tmp[g]; + } + } + + Tensor<1, dim *(dim + 1) / 2, Tensor<1, dim, Number>> result; + + // the diagonal part of Jacobian gradient comes first + for (unsigned int d = 0; d < dim; ++d) + for (unsigned int e = 0; e < dim; ++e) + result[d][e] = grad_jac_inv[d][d][e]; + + // then the upper-diagonal part + for (unsigned int d = 0, count = 0; d < dim; ++d) + for (unsigned int e = d + 1; e < dim; ++e, ++count) + for (unsigned int f = 0; f < dim; ++f) + result[dim + count][f] = grad_jac_inv[d][e][f]; + return result; + } + /** * Helper function called internally during the initialize function. */ @@ -670,9 +769,6 @@ namespace internal { FE_Nothing dummy_fe; - Tensor<3, dim, VectorizedArrayType> jac_grad, grad_jac_inv; - Tensor<1, dim, VectorizedArrayType> tmp; - // when we make comparisons about the size of Jacobians we need to // know the approximate size of typical entries in Jacobians. We need // to fix the Jacobian size once and for all. We choose the diameter @@ -713,33 +809,6 @@ namespace internal (update_flags & update_quadrature_points ? update_quadrature_points : update_default); - std::vector> n_q_points_1d(fe_values.size()), - step_size_cartesian(fe_values.size()); - for (unsigned int my_q = 0; my_q < fe_values.size(); ++my_q) - { - n_q_points_1d[my_q].resize( - mapping_info.cell_data[my_q].descriptor.size()); - step_size_cartesian[my_q].resize(n_q_points_1d[my_q].size()); - for (unsigned int hpq = 0; hpq < n_q_points_1d[my_q].size(); ++hpq) - { - n_q_points_1d[my_q][hpq] = mapping_info.cell_data[my_q] - .descriptor[hpq] - .quadrature_1d.size(); - - // To walk on the diagonal for lexicographic ordering, we have - // to jump one index ahead in each direction. For direction 0, - // this is just the next point, for direction 1, it means adding - // n_q_points_1d, and so on. - step_size_cartesian[my_q][hpq] = 0; - unsigned int factor = 1; - for (unsigned int d = 0; d < dim; ++d) - { - step_size_cartesian[my_q][hpq] += factor; - factor *= n_q_points_1d[my_q][hpq]; - } - } - } - const unsigned int end_cell = std::min(mapping_info.cell_type.size(), std::size_t(cell_range.second)); // loop over given cells @@ -876,110 +945,50 @@ namespace internal data.first[my_q].jacobians[0].push_back(inv_jac); if (update_flags & update_jacobian_grads) - { - // for second derivatives on the real cell, need - // also the gradient of the inverse Jacobian J. This - // involves some calculus and is done - // vectorized. This is very cheap compared to what - // fe_values does (in early 2011). If L is the - // gradient of the jacobian on the unit cell, the - // gradient of the inverse is given by - // (multidimensional calculus) - J * (J * L) * J - // (the third J is because we need to transform the - // gradient L from the unit to the real cell, and - // then apply the inverse Jacobian). Compare this - // with 1D with j(x) = 1/k(phi(x)), where j = phi' - // is the inverse of the jacobian and k is the - // derivative of the jacobian on the unit cell. Then - // j' = phi' k'/k^2 = j k' j^2. - - // compute: jac_grad = J*grad_unit(J^-1) - for (unsigned int d = 0; d < dim; ++d) - for (unsigned int e = 0; e < dim; ++e) - for (unsigned int f = 0; f < dim; ++f) - { - jac_grad[f][e][d] = - (inv_jac[f][0] * jacobian_grad[d][e][0]); - for (unsigned int g = 1; g < dim; ++g) - jac_grad[f][e][d] += - (inv_jac[f][g] * jacobian_grad[d][e][g]); - } - - // compute: transpose (-jac * jac_grad[d] * jac) - for (unsigned int d = 0; d < dim; ++d) - for (unsigned int e = 0; e < dim; ++e) - { - for (unsigned int f = 0; f < dim; ++f) - { - tmp[f] = VectorizedArrayType(); - for (unsigned int g = 0; g < dim; ++g) - tmp[f] -= - jac_grad[d][f][g] * inv_jac[g][e]; - } - - // needed for non-diagonal part of Jacobian - // grad - for (unsigned int f = 0; f < dim; ++f) - { - grad_jac_inv[f][d][e] = - inv_jac[f][0] * tmp[0]; - for (unsigned int g = 1; g < dim; ++g) - grad_jac_inv[f][d][e] += - inv_jac[f][g] * tmp[g]; - } - } - - // the diagonal part of Jacobian gradient comes first - Tensor<1, - dim *(dim + 1) / 2, - Tensor<1, dim, VectorizedArrayType>> - final_grad; - for (unsigned int d = 0; d < dim; ++d) - for (unsigned int e = 0; e < dim; ++e) - final_grad[d][e] = grad_jac_inv[d][d][e]; - - // then the upper-diagonal part - for (unsigned int d = 0, count = 0; d < dim; ++d) - for (unsigned int e = d + 1; e < dim; ++e, ++count) - for (unsigned int f = 0; f < dim; ++f) - final_grad[dim + count][f] = - grad_jac_inv[d][e][f]; - data.first[my_q].jacobian_gradients[0].push_back( - final_grad); - } + data.first[my_q].jacobian_gradients[0].push_back( + process_jacobian_gradient(inv_jac, jacobian_grad)); } } if (update_flags & update_quadrature_points) { // eventually we turn to the quadrature points that we can - // compress in case we have Cartesian cells. we also need to + // compress in case we have affine cells. we also need to // reorder them into arrays of vectorized data types. first // go through the cells and find out how much memory we need // to allocate for the quadrature points. We store - // n_q_points for all cells but Cartesian cells. For - // Cartesian cells, only need to store the values on a - // diagonal through the cell (n_q_points_1d). This will give - // (with some little indexing) the location of all - // quadrature points. + // n_q_points for general cells and a single value for + // Cartesian and affine cells (the position of the (0,0) + // point from the reference coordinates) const unsigned int old_size = data.first[my_q].quadrature_points.size(); data.first[my_q].quadrature_point_offsets.push_back(old_size); - if (mapping_info.get_cell_type(cell) == cartesian) + if (mapping_info.get_cell_type(cell) < general) { - for (unsigned int q = 0; - q < n_q_points_1d[my_q][fe_index]; - ++q) + Point quad_point; + for (unsigned int v = 0; + v < VectorizedArrayType::n_array_elements; + ++v) { - Point quad_point; + typename dealii::Triangulation::cell_iterator + cell_it( + &tria, + cells[cell * + VectorizedArrayType::n_array_elements + + v] + .first, + cells[cell * + VectorizedArrayType::n_array_elements + + v] + .second); + const Point p = + mapping.transform_unit_to_real_cell(cell_it, + Point()); for (unsigned int d = 0; d < dim; ++d) - quad_point[d] = - cell_data.quadrature_points - [q * step_size_cartesian[my_q][fe_index]][d]; - data.first[my_q].quadrature_points.push_back( - quad_point); + quad_point[d][v] = p[d]; } + data.first[my_q].quadrature_points.push_back(quad_point); } else { @@ -1038,8 +1047,10 @@ namespace internal { const unsigned int cell = lcell + first_cell; data_cells.data_index_offsets[cell] = - cell_type[cell] <= static_cast(affine) ? - indices_compressed[data_cells_local.data_index_offsets[lcell]] : + cell_type[cell] <= affine ? + (dim == structdim ? 2 : 1) * + indices_compressed[data_cells_local + .data_index_offsets[lcell]] : data_cells_local.data_index_offsets[lcell] + data_shift[0]; if (data_cells_local.quadrature_point_offsets.size() > lcell) data_cells.quadrature_point_offsets[cell] = @@ -1091,6 +1102,170 @@ namespace internal } } + + + /** + * This evaluates the mapping information on a range of cells calling + * into the tensor product interpolators of the matrix-free framework, + * using a polynomial expansion of the cell geometry in terms of + * MappingQ. + */ + template + void + compute_range_mapping_q( + const unsigned int begin_cell, + const unsigned int end_cell, + const std::vector & cell_type, + const std::vector & process_cell, + const UpdateFlags update_flags_cells, + const AlignedVector & plain_quadrature_points, + const ShapeInfo &shape_info, + MappingInfoStorage &my_data) + { + constexpr unsigned int n_lanes = VectorizedArrayType::n_array_elements; + constexpr unsigned int n_lanes_d = VectorizedDouble::n_array_elements; + + const unsigned int n_q_points = my_data.descriptor[0].n_q_points; + const unsigned int n_mapping_points = + shape_info.dofs_per_component_on_cell; + constexpr unsigned int hess_dim = dim * (dim + 1) / 2; + + AlignedVector cell_points(dim * n_mapping_points); + AlignedVector cell_quads(dim * n_q_points); + AlignedVector cell_grads(dim * dim * n_q_points); + AlignedVector cell_grad_grads(dim * hess_dim * + n_q_points); + AlignedVector scratch_data( + dim * (2 * n_q_points + 3 * n_mapping_points)); + + for (unsigned int cell = begin_cell; cell < end_cell; ++cell) + for (unsigned vv = 0; vv < n_lanes; vv += n_lanes_d) + { + if (cell_type[cell] > affine || process_cell[cell]) + { + unsigned int + start_indices[VectorizedDouble::n_array_elements]; + for (unsigned int v = 0; v < n_lanes_d; ++v) + start_indices[v] = + (cell * n_lanes + vv + v) * n_mapping_points * dim; + vectorized_load_and_transpose(n_mapping_points * dim, + plain_quadrature_points.data(), + start_indices, + cell_points.data()); + + SelectEvaluator::evaluate( + shape_info, + cell_points.data(), + cell_quads.data(), + cell_grads.data(), + cell_grad_grads.data(), + scratch_data.data(), + true, + true, + update_flags_cells & update_jacobian_grads); + } + if (update_flags_cells & update_quadrature_points) + { + Point *quadrature_points = + my_data.quadrature_points.data() + + my_data.quadrature_point_offsets[cell]; + if (cell_type[cell] <= affine) + for (unsigned int d = 0; d < dim; ++d) + for (unsigned int v = 0; v < n_lanes_d; ++v) + quadrature_points[0][d][vv + v] = + plain_quadrature_points + [(dim * (cell * n_lanes + vv + v) + d) * + n_mapping_points]; + else + for (unsigned int d = 0; d < dim; ++d) + for (unsigned int q = 0; q < n_q_points; ++q) + store_vectorized_array(cell_quads[q + d * n_q_points], + vv, + quadrature_points[q][d]); + } + + const unsigned int n_points = + cell_type[cell] <= affine ? 1 : n_q_points; + if (process_cell[cell]) + for (unsigned int q = 0; q < n_points; ++q) + { + const unsigned int idx = + my_data.data_index_offsets[cell] + q; + Tensor<2, dim, VectorizedDouble> jac; + for (unsigned int d = 0; d < dim; ++d) + for (unsigned int e = 0; e < dim; ++e) + jac[d][e] = cell_grads[q + (d * dim + e) * n_q_points]; + + // eliminate roundoff errors + if (cell_type[cell] == cartesian) + for (unsigned int d = 0; d < dim; ++d) + for (unsigned int e = 0; e < dim; ++e) + if (d != e) + jac[d][e] = 0.; + + const VectorizedDouble jac_det = determinant(jac); + const Tensor<2, dim, VectorizedDouble> inv_jac = + transpose(invert(jac)); + + if (cell_type[cell] <= affine) + { + store_vectorized_array(jac_det, + vv, + my_data.JxW_values[idx]); + + for (unsigned int d = 0; d < dim; ++d) + for (unsigned int e = 0; e < dim; ++e) + store_vectorized_array( + jac[d][e], + vv, + my_data.jacobians[0][idx + 1][d][e]); + } + else + { + const double weight = + my_data.descriptor[0].quadrature.weight(q); + store_vectorized_array(jac_det * weight, + vv, + my_data.JxW_values[idx]); + } + for (unsigned int d = 0; d < dim; ++d) + for (unsigned int e = 0; e < dim; ++e) + store_vectorized_array(inv_jac[d][e], + vv, + my_data.jacobians[0][idx][d][e]); + + if (update_flags_cells & update_jacobian_grads && + cell_type[cell] > affine) + { + Tensor<3, dim, VectorizedDouble> jac_grad; + for (unsigned int d = 0; d < dim; ++d) + { + for (unsigned int e = 0; e < dim; ++e) + jac_grad[d][e][e] = + cell_grad_grads[q + (d * hess_dim + e) * + n_q_points]; + for (unsigned int c = dim, e = 0; e < dim; ++e) + for (unsigned int f = e + 1; f < dim; ++f, ++c) + jac_grad[d][e][f] = jac_grad[d][f][e] = + cell_grad_grads[q + (d * hess_dim + c) * + n_q_points]; + const auto inv_jac_grad = + process_jacobian_gradient(inv_jac, jac_grad); + for (unsigned int d = 0; d < hess_dim; ++d) + for (unsigned int e = 0; e < dim; ++e) + store_vectorized_array( + inv_jac_grad[d][e], + vv, + my_data.jacobian_gradients[0][idx][d][e]); + } + } + } + } + } + } // namespace ExtractCellHelper @@ -1174,7 +1349,7 @@ namespace internal cell_data[my_q].data_index_offsets.resize(cell_type.size()); std::vector> shift( data_cells_local.size()); - shift[0][0] = n_constant_jacobians; + shift[0][0] = 2 * n_constant_jacobians; shift[0][1] = 0; for (unsigned int i = 1; i < data_cells_local.size(); ++i) { @@ -1215,7 +1390,11 @@ namespace internal cell_data[my_q]); // finally, insert the constant cell data at the beginning (the - // other tasks can already start copying the non-constant data) + // other tasks can already start copying the non-constant + // data). Note that we use two slots for the constant data to + // accommodate for both the inverse transposed Jacobian (that we + // need for derivatives) and the Jacobian (that we need for + // quadrature points) if (my_q == 0) { for (const auto &it : data_cells_local[0].second.data) @@ -1227,16 +1406,18 @@ namespace internal ++v) jac[d][e][v] = it.first[d][e][v]; AssertIndexRange(it.second, n_constant_jacobians); - const std::size_t index = it.second; - cell_data[my_q].JxW_values[index] = determinant(jac); + const std::size_t index = it.second; + cell_data[my_q].JxW_values[2 * index] = determinant(jac); // invert and transpose jac - cell_data[my_q].jacobians[0][index] = transpose(invert(jac)); + cell_data[my_q].jacobians[0][2 * index] = + transpose(invert(jac)); + cell_data[my_q].jacobians[0][2 * index + 1] = jac; // second derivative of transformation is zero on affine cells } } else { - for (unsigned int i = 0; i < n_constant_jacobians; ++i) + for (unsigned int i = 0; i < 2 * n_constant_jacobians; ++i) { cell_data[my_q].JxW_values[i] = cell_data[0].JxW_values[i]; cell_data[my_q].jacobians[0][i] = @@ -1738,6 +1919,238 @@ namespace internal } } + + + /** + * This evaluates the mapping information on a range of cells calling + * into the tensor product interpolators of the matrix-free framework, + * using a polynomial expansion of the cell geometry in terms of + * MappingQ. + */ + template + void + compute_range_mapping_q( + const unsigned int begin_face, + const unsigned int end_face, + const std::vector< + FaceToCellTopology> &faces, + const std::vector & face_type, + const std::vector & process_face, + const UpdateFlags update_flags_faces, + const AlignedVector & plain_quadrature_points, + const ShapeInfo &shape_info, + MappingInfoStorage &my_data) + { + constexpr unsigned int n_lanes = VectorizedArrayType::n_array_elements; + constexpr unsigned int n_lanes_d = VectorizedDouble::n_array_elements; + + const unsigned int n_q_points = my_data.descriptor[0].n_q_points; + const unsigned int n_mapping_points = + shape_info.dofs_per_component_on_cell; + + AlignedVector cell_points(dim * n_mapping_points); + AlignedVector face_quads(dim * n_q_points); + AlignedVector face_grads(dim * dim * n_q_points); + AlignedVector scratch_data( + dim * (2 * n_q_points + 3 * n_mapping_points)); + + for (unsigned int face = begin_face; face < end_face; ++face) + for (unsigned vv = 0; vv < n_lanes; vv += n_lanes_d) + { + // load the geometry field for all SIMD lanes + unsigned int start_indices[VectorizedDouble::n_array_elements]; + const unsigned int face_no = faces[face].interior_face_no; + for (unsigned int v = 0; v < n_lanes_d; ++v) + if (faces[face].cells_interior[vv + v] != + numbers::invalid_unsigned_int) + start_indices[v] = + faces[face].cells_interior[vv + v] * n_mapping_points * dim; + else + start_indices[v] = + faces[face].cells_interior[0] * n_mapping_points * dim; + vectorized_load_and_transpose(n_mapping_points * dim, + plain_quadrature_points.data(), + start_indices, + cell_points.data()); + + // now let the matrix-free evaluators provide us with the + // data on faces + FEFaceEvaluationSelector:: + evaluate(shape_info, + cell_points.data(), + face_quads.data(), + face_grads.data(), + scratch_data.data(), + true, + true, + face_no, + GeometryInfo::max_children_per_cell, + faces[face].face_orientation > 8 ? + faces[face].face_orientation - 8 : + 0, + my_data.descriptor[0].face_orientations); + + + if (update_flags_faces & update_quadrature_points) + for (unsigned int q = 0; q < n_q_points; ++q) + for (unsigned int d = 0; d < dim; ++d) + store_vectorized_array( + face_quads[d * n_q_points + q], + vv, + my_data.quadrature_points + [my_data.quadrature_point_offsets[face] + q][d]); + + if (process_face[face] == false) + continue; + + // go through the faces and fill the result + const unsigned int offset = my_data.data_index_offsets[face]; + const unsigned int n_points_compute = + face_type[face] <= affine ? 1 : n_q_points; + for (unsigned int q = 0; q < n_points_compute; ++q) + { + Tensor<2, dim, VectorizedDouble> jac; + for (unsigned int e = 0; e < dim; ++e) + { + const unsigned int ee = + ExtractFaceHelper::reorder_face_derivative_indices( + face_no, e); + for (unsigned int d = 0; d < dim; ++d) + jac[d][ee] = face_grads[(d * dim + e) * n_q_points + q]; + } + Tensor<2, dim, VectorizedDouble> inv_jac = invert(jac); + for (unsigned int e = 0; e < dim; ++e) + { + const unsigned int ee = + ExtractFaceHelper::reorder_face_derivative_indices( + face_no, e); + for (unsigned int d = 0; d < dim; ++d) + store_vectorized_array( + inv_jac[ee][d], + vv, + my_data.jacobians[0][offset + q][d][e]); + } + + std::array, dim - 1> + tangential_vectors; + for (unsigned int d = 0; d != dim - 1; ++d) + for (unsigned int e = 0; e < dim; ++e) + for (unsigned int f = 0; f < dim; ++f) + tangential_vectors[d][e] += + jac[e][f] * + GeometryInfo::unit_tangential_vectors[face_no][d] + [f]; + + Tensor<1, dim, VectorizedDouble> boundary_form; + if (dim == 1) + boundary_form[0] = face_no == 0 ? -1. : 1.; + else if (dim == 2) + boundary_form = cross_product_2d(tangential_vectors[0]); + else if (dim == 3) + boundary_form = cross_product_3d(tangential_vectors[0], + tangential_vectors[1]); + else + Assert(false, ExcNotImplemented()); + + const VectorizedDouble JxW = + boundary_form.norm() * + (face_type[face] <= affine ? + 1. : + my_data.descriptor[0].quadrature.weight(q)); + + store_vectorized_array(JxW, + vv, + my_data.JxW_values[offset + q]); + + const Tensor<1, dim, VectorizedDouble> normal = + boundary_form / boundary_form.norm(); + + for (unsigned int d = 0; d < dim; ++d) + store_vectorized_array( + normal[d], vv, my_data.normal_vectors[offset + q][d]); + + my_data.normals_times_jacobians[0][offset + q] = + my_data.normal_vectors[offset + q] * + my_data.jacobians[0][offset + q]; + } + + if (faces[face].cells_exterior[0] != + numbers::invalid_unsigned_int) + { + for (unsigned int v = 0; v < n_lanes_d; ++v) + if (faces[face].cells_exterior[vv + v] != + numbers::invalid_unsigned_int) + start_indices[v] = faces[face].cells_exterior[vv + v] * + n_mapping_points * dim; + else + start_indices[v] = + faces[face].cells_exterior[0] * n_mapping_points * dim; + + vectorized_load_and_transpose(n_mapping_points * dim, + plain_quadrature_points.data(), + start_indices, + cell_points.data()); + + FEFaceEvaluationSelector:: + evaluate(shape_info, + cell_points.data(), + face_quads.data(), + face_grads.data(), + scratch_data.data(), + false, + true, + faces[face].exterior_face_no, + faces[face].subface_index, + faces[face].face_orientation < 8 ? + faces[face].face_orientation : + 0, + my_data.descriptor[0].face_orientations); + + for (unsigned int q = 0; q < n_points_compute; ++q) + { + Tensor<2, dim, VectorizedDouble> jac; + for (unsigned int e = 0; e < dim; ++e) + { + const unsigned int ee = + ExtractFaceHelper::reorder_face_derivative_indices< + dim>(faces[face].exterior_face_no, e); + for (unsigned int d = 0; d < dim; ++d) + jac[d][ee] = + face_grads[(d * dim + e) * n_q_points + q]; + } + Tensor<2, dim, VectorizedDouble> inv_jac = invert(jac); + for (unsigned int e = 0; e < dim; ++e) + { + const unsigned int ee = + ExtractFaceHelper::reorder_face_derivative_indices< + dim>(faces[face].exterior_face_no, e); + for (unsigned int d = 0; d < dim; ++d) + store_vectorized_array( + inv_jac[ee][d], + vv, + my_data.jacobians[1][offset + q][d][e]); + } + my_data.normals_times_jacobians[1][offset + q] = + my_data.normal_vectors[offset + q] * + my_data.jacobians[1][offset + q]; + } + } + } + } + } // namespace ExtractFaceHelper @@ -1755,7 +2168,8 @@ namespace internal if (faces.size() == 0) return; - // Create as many chunks of cells as we have threads and spawn the work + // Create as many chunks of cells as we have threads and spawn the + // work unsigned int work_per_chunk = std::max(std::size_t(8), (faces.size() + MultithreadInfo::n_threads() - 1) / @@ -1767,8 +2181,8 @@ namespace internal ExtractFaceHelper:: CompressedFaceData>> data_faces_local; - // Reserve enough space to avoid re-allocation (which would destroy the - // references passed to the tasks!) + // Reserve enough space to avoid re-allocation (which would destroy + // the references passed to the tasks!) data_faces_local.reserve(MultithreadInfo::n_threads()); { @@ -1937,6 +2351,414 @@ namespace internal + template + void + MappingInfo::compute_mapping_q( + const dealii::Triangulation & tria, + const std::vector> &cell_array, + const std::vector< + FaceToCellTopology> &faces) + { + // step 1: extract quadrature point data with the data appropriate for + // MappingQGeneric + const MappingQGeneric *mapping_q = + dynamic_cast *>(&*this->mapping); + Assert(mapping_q != nullptr, ExcInternalError()); + + const unsigned int mapping_degree = mapping_q->get_degree(); + const unsigned int n_mapping_points = + Utilities::pow(mapping_degree + 1, dim); + AlignedVector plain_quadrature_points(cell_array.size() * + n_mapping_points * dim); + + const double jacobian_size = ExtractCellHelper::get_jacobian_size(tria); + + std::vector cell_data_index(cell_array.size()); + std::vector preliminary_cell_type(cell_array.size()); + { + FE_Nothing dummy_fe; + QGaussLobatto quadrature(mapping_degree + 1); + + FEValues fe_values(*mapping_q, + dummy_fe, + quadrature, + update_quadrature_points | update_jacobians); + + // we include a map to store some compressed information about the + // Jacobians which we collect by a stencil-like pattern around the + // first quadrature point on the cell - we use a relatively coarse + // tolerance to account for some inaccuracies in the manifold + // evaluation + const FPArrayComparator comparator(1e4 * jacobian_size); + std::map, dim + 1>, + unsigned int, + FPArrayComparator> + compressed_jacobians(comparator); + + unsigned int n_data_buckets = 0; + for (unsigned int cell = 0; cell < cell_array.size(); ++cell) + { + typename dealii::Triangulation::cell_iterator cell_it( + &tria, cell_array[cell].first, cell_array[cell].second); + fe_values.reinit(cell_it); + for (unsigned int d = 0; d < dim; ++d) + for (unsigned int q = 0; q < n_mapping_points; ++q) + plain_quadrature_points[(cell * dim + d) * n_mapping_points + + q] = fe_values.quadrature_point(q)[d]; + + // store the first, second, n-th and n^2-th one along a + // stencil-like pattern + std::array, dim + 1> jacobians_on_stencil; + jacobians_on_stencil[0] = + Tensor<2, dim, double>(fe_values.jacobian(0)); + for (unsigned int d = 0, skip = 1; d < dim; + ++d, skip *= (mapping_degree + 1)) + jacobians_on_stencil[1 + d] = + Tensor<2, dim, double>(fe_values.jacobian(skip)); + + // check in the map for the index of this cell + auto inserted = compressed_jacobians.insert( + std::make_pair(jacobians_on_stencil, cell)); + bool add_this_cell = inserted.second; + if (inserted.second == false) + { + // check if the found duplicate really is a translation and + // the similarity identified by the map is not by accident + double max_distance = 0; + const double *ptr_origin = + plain_quadrature_points.data() + + inserted.first->second * dim * n_mapping_points; + const double *ptr_mine = plain_quadrature_points.data() + + cell * dim * n_mapping_points; + for (unsigned int d = 0; d < dim; ++d) + { + const double translate_d = + ptr_origin[d * n_mapping_points] - + ptr_mine[d * n_mapping_points]; + for (unsigned int q = 1; q < n_mapping_points; ++q) + max_distance = + std::max(std::abs(ptr_origin[d * n_mapping_points + q] - + ptr_mine[d * n_mapping_points + q] - + translate_d), + max_distance); + } + + // this is not a duplicate, must add it again + if (max_distance > 1e-10 * jacobian_size) + add_this_cell = true; + } + + if (add_this_cell == true) + { + // check whether cell is Cartesian/affine/general + GeometryType type = cartesian; + for (unsigned int d = 0; d < dim; ++d) + for (unsigned int e = 0; e < dim; ++e) + if (d != e) + if (std::abs(inserted.first->first[0][d][e]) > + 1e-12 * jacobian_size) + type = affine; + + for (unsigned int q = 1; q < n_mapping_points; ++q) + for (unsigned int d = 0; d < dim; ++d) + for (unsigned int e = 0; e < dim; ++e) + if (std::abs(fe_values.jacobian(q)[d][e] - + fe_values.jacobian(0)[d][e]) > + 1e-12 * jacobian_size) + { + type = general; + goto endloop; + } + endloop: + cell_data_index[cell] = n_data_buckets; + preliminary_cell_type[cell] = type; + ++n_data_buckets; + } + else + { + cell_data_index[cell] = cell_data_index[inserted.first->second]; + preliminary_cell_type[cell] = + preliminary_cell_type[inserted.first->second]; + } + } + } + + // step 2: compute the appropriate evaluation matrices for cells and + // faces + + // We want to use vectorization for computing the quantities, but must + // evaluate the geometry in double precision; thus, for floats we need + // to do things in two sweeps and convert the final result. + constexpr unsigned int n_lanes = VectorizedArrayType::n_array_elements; + using VectorizedDouble = + VectorizedArray::value && n_lanes > 1) ? + n_lanes / 2 : + n_lanes)>; + + // Create a ShapeInfo object to provide the necessary interpolators to + // the various quadrature points. Note that it is initialized with the + // finite element fe_geometry using the degree of the mapping, which is + // not the same as the degree of the underlying finite element shape + // functions or the quadrature points; shape info is merely a vehicle to + // return us the right interpolation matrices from the cell support + // points to the cell and face quadrature points. + std::vector> shape_infos(cell_data.size()); + { + FE_DGQ fe_geometry(mapping_degree); + for (unsigned int my_q = 0; my_q < cell_data.size(); ++my_q) + shape_infos[my_q].reinit(cell_data[my_q].descriptor[0].quadrature_1d, + fe_geometry); + } + + // step 3: find compression of cells with vectorization + std::map, unsigned int> compressed_data; + + cell_type.resize(cell_array.size() / n_lanes); + std::vector process_cell(cell_type.size()); + std::vector cell_data_index_vect(cell_type.size()); + + for (unsigned int cell = 0; cell < cell_array.size(); cell += n_lanes) + { + std::pair, unsigned int> + data_indices; + for (unsigned int i = 0; i < n_lanes; ++i) + data_indices.first[i] = cell_data_index[cell + i]; + data_indices.second = cell / n_lanes; + + auto inserted = compressed_data.insert(data_indices); + + process_cell[cell / n_lanes] = inserted.second; + if (inserted.second == true) + cell_data_index_vect[cell / n_lanes] = data_indices.second; + else + cell_data_index_vect[cell / n_lanes] = inserted.first->second; + + cell_type[cell / n_lanes] = + *std::max_element(preliminary_cell_type.data() + cell, + preliminary_cell_type.data() + cell + n_lanes); + } + + // step 4: compute the data on cells from the cached quadrature + // points, filling up all SIMD lanes as appropriate + for (unsigned int my_q = 0; my_q < cell_data.size(); ++my_q) + { + MappingInfoStorage &my_data = + cell_data[my_q]; + + // step 4a: set the index offsets, find out how much to allocate, + // and allocate the memory + const unsigned int n_q_points = my_data.descriptor[0].n_q_points; + unsigned int max_size = 0; + my_data.data_index_offsets.resize(cell_type.size()); + for (unsigned int cell = 0; cell < cell_type.size(); ++cell) + { + if (process_cell[cell] == false) + my_data.data_index_offsets[cell] = + my_data.data_index_offsets[cell_data_index_vect[cell]]; + else + my_data.data_index_offsets[cell] = max_size; + max_size = + std::max(max_size, + my_data.data_index_offsets[cell] + + (cell_type[cell] <= affine ? 2 : n_q_points)); + } + + my_data.JxW_values.resize_fast(max_size); + my_data.jacobians[0].resize_fast(max_size); + if (update_flags_cells & update_jacobian_grads) + my_data.jacobian_gradients[0].resize_fast(max_size); + + if (update_flags_cells & update_quadrature_points) + { + my_data.quadrature_point_offsets.resize(cell_type.size()); + for (unsigned int cell = 1; cell < cell_type.size(); ++cell) + if (cell_type[cell - 1] <= affine) + my_data.quadrature_point_offsets[cell] = + my_data.quadrature_point_offsets[cell - 1] + 1; + else + my_data.quadrature_point_offsets[cell] = + my_data.quadrature_point_offsets[cell - 1] + n_q_points; + my_data.quadrature_points.resize_fast( + my_data.quadrature_point_offsets.back() + + (cell_type.back() <= affine ? 1 : n_q_points)); + } + + // step 4b: go through the cells and compute the information using + // similar evaluators as for the matrix-free integrals + ExtractCellHelper::compute_range_mapping_q( + 0, + cell_type.size(), + cell_type, + process_cell, + update_flags_cells, + plain_quadrature_points, + shape_infos[my_q], + my_data); + } + + if (faces.empty()) + return; + + // step 5: find compression of faces with vectorization + std::map, unsigned int> + compressed_faces; + + face_type.resize(faces.size()); + std::vector process_face(face_type.size()); + std::vector face_data_index_vect(face_type.size()); + + for (unsigned int face = 0; face < faces.size(); ++face) + { + std::pair, unsigned int> + data_indices; + for (unsigned int i = 0; i < n_lanes; ++i) + if (faces[face].cells_interior[i] != numbers::invalid_unsigned_int) + data_indices.first[i] = + cell_data_index[faces[face].cells_interior[i]]; + else + data_indices.first[i] = data_indices.first[0]; + for (unsigned int i = 0; i < n_lanes; ++i) + data_indices.first[n_lanes + i] = data_indices.first[i]; + for (unsigned int i = 0; i < n_lanes; ++i) + if (faces[face].cells_exterior[i] != numbers::invalid_unsigned_int) + data_indices.first[n_lanes + i] = + cell_data_index[faces[face].cells_exterior[i]]; + data_indices.first[2 * n_lanes] = faces[face].interior_face_no; + data_indices.first[2 * n_lanes + 1] = faces[face].exterior_face_no; + data_indices.first[2 * n_lanes + 2] = faces[face].subface_index; + + data_indices.second = face; + + auto inserted = compressed_faces.insert(data_indices); + + process_face[face] = inserted.second; + if (inserted.second == true) + face_data_index_vect[face] = face; + else + face_data_index_vect[face] = inserted.first->second; + + face_type[face] = cartesian; + for (unsigned int i = 0; i < n_lanes; ++i) + if (faces[face].cells_interior[i] != numbers::invalid_unsigned_int) + face_type[face] = + std::max(face_type[face], + preliminary_cell_type[faces[face].cells_interior[i]]); + for (unsigned int i = 0; i < n_lanes; ++i) + if (faces[face].cells_exterior[i] != numbers::invalid_unsigned_int) + face_type[face] = + std::max(face_type[face], + preliminary_cell_type[faces[face].cells_exterior[i]]); + } + + // step 6: compute the data on faces from the cached cell quadrature + // points, filling up all SIMD lanes as appropriate + for (unsigned int my_q = 0; my_q < face_data.size(); ++my_q) + { + MappingInfoStorage + &my_data = face_data[my_q]; + + // step 6a: set the index offsets, find out how much to allocate, + // and allocate the memory + const unsigned int n_q_points = my_data.descriptor[0].n_q_points; + unsigned int max_size = 0; + my_data.data_index_offsets.resize(face_type.size()); + for (unsigned int face = 0; face < face_type.size(); ++face) + { + if (process_face[face] == false) + my_data.data_index_offsets[face] = + my_data.data_index_offsets[face_data_index_vect[face]]; + else + my_data.data_index_offsets[face] = max_size; + max_size = + std::max(max_size, + my_data.data_index_offsets[face] + + (face_type[face] <= affine ? 1 : n_q_points)); + } + + const UpdateFlags update_flags_common = + update_flags_boundary_faces | update_flags_inner_faces; + + my_data.JxW_values.resize_fast(max_size); + my_data.normal_vectors.resize_fast(max_size); + my_data.jacobians[0].resize_fast(max_size); + my_data.jacobians[1].resize_fast(max_size); + if (update_flags_common & update_jacobian_grads) + { + my_data.jacobian_gradients[0].resize_fast(max_size); + my_data.jacobian_gradients[1].resize_fast(max_size); + } + my_data.normals_times_jacobians[0].resize_fast(max_size); + my_data.normals_times_jacobians[1].resize_fast(max_size); + + if (update_flags_cells & update_quadrature_points) + { + my_data.quadrature_point_offsets.resize(face_type.size()); + my_data.quadrature_point_offsets[0] = 0; + for (unsigned int face = 1; face < faces.size(); ++face) + my_data.quadrature_point_offsets[face] = + n_q_points + my_data.quadrature_point_offsets[face - 1]; + my_data.quadrature_points.resize_fast(face_type.size() * + n_q_points); + } + + // step 6b: go through the faces and compute the information using + // similar evaluators as for the matrix-free face integrals + ExtractFaceHelper::compute_range_mapping_q( + 0, + face_type.size(), + faces, + face_type, + process_face, + update_flags_common, + plain_quadrature_points, + shape_infos[my_q], + my_data); + } + + // step 6c: figure out if normal vectors are the same on some of the + // faces which allows us to set the flat_faces face type + unsigned int quad_with_most_points = 0; + for (unsigned int my_q = 1; my_q < face_data.size(); ++my_q) + if (face_data[my_q].descriptor[0].n_q_points > + face_data[quad_with_most_points].descriptor[0].n_q_points) + quad_with_most_points = my_q; + for (unsigned int face = 0; face < face_type.size(); ++face) + if (face_type[face] == general) + { + const unsigned int n_q_points = + face_data[quad_with_most_points].descriptor[0].n_q_points; + const Tensor<1, dim, VectorizedArrayType> *normals = + face_data[quad_with_most_points].normal_vectors.data() + + face_data[quad_with_most_points].data_index_offsets[face]; + VectorizedArrayType distance = 0.; + for (unsigned int q = 1; q < n_q_points; ++q) + distance += (normals[q] - normals[0]).norm_square(); + bool all_small = true; + for (unsigned int v = 0; v < n_lanes; ++v) + if (distance[v] > 50. * std::numeric_limits::epsilon() * + std::numeric_limits::epsilon() * + n_q_points) + all_small = false; + if (all_small) + face_type[face] = flat_faces; + } + + // step 7: compute the face data by cells. This still needs to be + // transitioned to extracting the information from cell quadrature + // points but we need to figure out the correct indices of neighbors + // within the list of arrays still + initialize_faces_by_cells(tria, cell_array, *this->mapping); + } + + + template void MappingInfo::initialize_faces_by_cells( @@ -2327,7 +3149,26 @@ namespace internal return false; } - } // end of namespace MatrixFreeFunctions + + + template + template + bool + FPArrayComparator:: + operator()(const std::array, dim + 1> &t1, + const std::array, dim + 1> &t2) const + { + for (unsigned int i = 0; i < t1.size(); ++i) + for (unsigned int d = 0; d < dim; ++d) + for (unsigned int e = 0; e < dim; ++e) + if (t1[i][d][e] < t2[i][d][e] - tolerance) + return true; + else if (t1[i][d][e] > t2[i][d][e] + tolerance) + return false; + return false; + } + + } // namespace MatrixFreeFunctions } // end of namespace internal DEAL_II_NAMESPACE_CLOSE diff --git a/tests/matrix_free/compress_mapping.cc b/tests/matrix_free/compress_mapping.cc index 8493156b99..cc9cf159cf 100644 --- a/tests/matrix_free/compress_mapping.cc +++ b/tests/matrix_free/compress_mapping.cc @@ -123,10 +123,9 @@ test_cube() for (unsigned int i = 0; i < n_macro_cells; ++i) n_cell_types[mf.get_mapping_info().get_cell_type(i)]++; - // should have one Cartesian cell and no other - // cell type + // should have one Cartesian cell and no other cell type AssertDimension(n_cell_types[0], n_macro_cells); - AssertDimension(mf.get_mapping_info().cell_data[0].jacobians[0].size(), 1); + AssertDimension(mf.get_mapping_info().cell_data[0].jacobians[0].size(), 2); Assert(n_macro_cells > 1, ExcInternalError()); deallog << "OK" << std::endl; } @@ -170,7 +169,7 @@ test_parallelogram() // should have one affine cell and no other // cell type AssertDimension(n_cell_types[1], n_macro_cells); - AssertDimension(mf.get_mapping_info().cell_data[0].jacobians[0].size(), 1); + AssertDimension(mf.get_mapping_info().cell_data[0].jacobians[0].size(), 2); Assert(n_macro_cells > 1, ExcInternalError()); deallog << "OK" << std::endl; }