From: Martin Kronbichler Date: Sat, 7 Apr 2018 21:15:28 +0000 (+0200) Subject: Convert data structures in MF::MappingInfo to a face-compatible variant. X-Git-Tag: v9.0.0-rc1~194^2~1 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=11584e6e0aaf1aa435f9cd2751ee252ff4647f4a;p=dealii.git Convert data structures in MF::MappingInfo to a face-compatible variant. --- diff --git a/include/deal.II/matrix_free/dof_info.templates.h b/include/deal.II/matrix_free/dof_info.templates.h index 0f0062c3e5..894167a1ca 100644 --- a/include/deal.II/matrix_free/dof_info.templates.h +++ b/include/deal.II/matrix_free/dof_info.templates.h @@ -23,6 +23,7 @@ #include #include #include +#include DEAL_II_NAMESPACE_OPEN diff --git a/include/deal.II/matrix_free/fe_evaluation.h b/include/deal.II/matrix_free/fe_evaluation.h index d4e83961fb..8b6a6579fa 100644 --- a/include/deal.II/matrix_free/fe_evaluation.h +++ b/include/deal.II/matrix_free/fe_evaluation.h @@ -135,17 +135,23 @@ public: void reinit (const typename Triangulation::cell_iterator &cell); /** - * For the transformation information stored in MappingInfo, this function - * returns the index which belongs to the current cell as specified in @p - * reinit. Note that MappingInfo has different fields for Cartesian cells, - * cells with affine mapping and with general mappings, so in order to - * access the correct data, this interface must be used together with - * get_cell_type. + * @deprecated Use get_mapping_data_index_offset() instead. */ + DEAL_II_DEPRECATED unsigned int get_cell_data_number() const; /** - * Return the type of the cell the @p reinit function has been called for. + * Return the index offset within the geometry fields for the cell the @p + * reinit() function has been called for. This index can be used to access + * an index into a field that has the same compression behavior as the + * Jacobian of the geometry, e.g., to store an effective coefficient tensors + * that combines a coefficient with the geometry for lower memory transfer + * as the available data fields. + */ + unsigned int get_mapping_data_index_offset() const; + + /** + * Return the type of the cell the @p reinit() function has been called for. * Valid values are @p cartesian for Cartesian cells (which allows for * considerable data compression), @p affine for cells with affine mappings, * and @p general for general cells without any compressed storage applied. @@ -481,6 +487,15 @@ public: */ VectorizedArray JxW(const unsigned int q_point) const; + /** + * Gets the inverse and transposed version of Jacobian of the mapping + * between the unit to the real cell (representing the covariant + * transformation). This is exactly the matrix used internally to transform + * the unit cell gradients to gradients on the real cell. + */ + Tensor<2,dim,VectorizedArray > + inverse_jacobian(const unsigned int q_index) const; + //@} /** @@ -800,12 +815,12 @@ protected: const internal::MatrixFreeFunctions::DoFInfo *dof_info; /** - * Stores a pointer to the underlying transformation data from unit to real - * cells for the given quadrature formula specified at construction. Also - * contained in matrix_info, but it simplifies code if we store a reference - * to it. + * Stores a pointer to the underlying transformation data from unit to + * real cells for the given quadrature formula specified at construction. + * Also contained in matrix_info, but it simplifies code if we store a + * reference to it. */ - const internal::MatrixFreeFunctions::MappingInfo *mapping_info; + const internal::MatrixFreeFunctions::MappingInfoStorage *mapping_data; /** * Stores a pointer to the unit cell shape data, i.e., values, gradients and @@ -815,12 +830,6 @@ protected: */ const internal::MatrixFreeFunctions::ShapeInfo> *data; - /** - * A pointer to the Cartesian Jacobian information of the present cell. Only - * set to a useful value if on a Cartesian cell, otherwise zero. - */ - const Tensor<1,dim,VectorizedArray > *cartesian_data; - /** * A pointer to the Jacobian information of the present cell. Only set to a * useful value if on a non-Cartesian cell. @@ -838,26 +847,7 @@ protected: /** * A pointer to the quadrature weights of the underlying quadrature formula. */ - const VectorizedArray *quadrature_weights; - - /** - * A pointer to the quadrature points on the present cell. - */ - const Point > *quadrature_points; - - /** - * A pointer to the diagonal part of the Jacobian gradient on the present - * cell. Only set to a useful value if on a general cell with non-constant - * Jacobian. - */ - const Tensor<2,dim,VectorizedArray > *jacobian_grad; - - /** - * A pointer to the upper diagonal part of the Jacobian gradient on the - * present cell. Only set to a useful value if on a general cell with non- - * constant Jacobian. - */ - const Tensor<1,(dim>1?dim*(dim-1)/2:1),Tensor<1,dim,VectorizedArray > > * jacobian_grad_upper; + const Number *quadrature_weights; /** * After a call to reinit(), stores the number of the cell we are currently @@ -873,11 +863,6 @@ protected: */ internal::MatrixFreeFunctions::CellType cell_type; - /** - * The stride to access the correct data in MappingInfo. - */ - unsigned int cell_data_number; - /** * Debug information to track whether dof values have been initialized * before accessed. Used to control exceptions when uninitialized data is @@ -926,12 +911,6 @@ protected: */ std::shared_ptr > mapped_geometry; - /** - * For use with on-the-fly evaluation, provide a data structure to store the - * global dof indices on the current cell from a reinit call. - */ - std::vector old_style_dof_indices; - /** * For a FiniteElement with more than one finite element, select at which * component this data structure should start. @@ -2112,28 +2091,22 @@ FEEvaluationBase : 0), active_quad_index (fe_degree != numbers::invalid_unsigned_int ? - data_in.get_mapping_info(). - mapping_data_gen[quad_no_in]. - quad_index_from_n_q_points(n_q_points) + (data_in.get_mapping_info().cell_data[quad_no_in]. + quad_index_from_n_q_points(n_q_points)) : 0), matrix_info (&data_in), dof_info (&data_in.get_dof_info(fe_no_in)), - mapping_info (&data_in.get_mapping_info()), + mapping_data (internal::MatrixFreeFunctions::MappingInfoCellsOrFaces::get + (data_in.get_mapping_info(), quad_no)), data (&data_in.get_shape_info (fe_no_in, quad_no_in, active_fe_index, active_quad_index)), - cartesian_data (nullptr), jacobian (nullptr), J_value (nullptr), - quadrature_weights (mapping_info->mapping_data_gen[quad_no]. - quadrature_weights[active_quad_index].begin()), - quadrature_points (nullptr), - jacobian_grad (nullptr), - jacobian_grad_upper(nullptr), + quadrature_weights (mapping_data->descriptor[active_quad_index].quadrature_weights.begin()), cell (numbers::invalid_unsigned_int), - cell_type (internal::MatrixFreeFunctions::undefined), - cell_data_number (numbers::invalid_unsigned_int), + cell_type (internal::MatrixFreeFunctions::general), dof_values_initialized (false), values_quad_initialized (false), gradients_quad_initialized(false), @@ -2150,7 +2123,7 @@ FEEvaluationBase AssertDimension (data->dofs_per_component_on_cell*n_fe_components, dof_info->dofs_per_cell[active_fe_index]); AssertDimension (data->n_q_points, - mapping_info->mapping_data_gen[quad_no].n_q_points[active_quad_index]); + mapping_data->descriptor[active_quad_index].n_q_points); Assert (n_fe_components == 1 || n_components == 1 || n_components == n_fe_components, @@ -2183,19 +2156,14 @@ FEEvaluationBase active_quad_index (numbers::invalid_unsigned_int), matrix_info (nullptr), dof_info (nullptr), - mapping_info (nullptr), + mapping_data (nullptr), // select the correct base element from the given FE component data (new internal::MatrixFreeFunctions::ShapeInfo>(quadrature, fe, fe.component_to_base_index(first_selected_component).first)), - cartesian_data (nullptr), jacobian (nullptr), J_value (nullptr), quadrature_weights (nullptr), - quadrature_points (nullptr), - jacobian_grad (nullptr), - jacobian_grad_upper(nullptr), cell (0), cell_type (internal::MatrixFreeFunctions::general), - cell_data_number (numbers::invalid_unsigned_int), dof_values_initialized (false), values_quad_initialized (false), gradients_quad_initialized(false), @@ -2218,9 +2186,10 @@ FEEvaluationBase mapped_geometry = std::make_shared > (mapping, quadrature, update_flags); - jacobian = mapped_geometry->get_inverse_jacobians().begin(); - J_value = mapped_geometry->get_JxW_values().begin(); - quadrature_points = mapped_geometry->get_quadrature_points().begin(); + + mapping_data = &mapped_geometry->get_data_storage(); + jacobian = mapped_geometry->get_data_storage().jacobians[0].begin(); + J_value = mapped_geometry->get_data_storage().JxW_values.begin(); Assert(fe.element_multiplicity(base_element_number) == 1 || fe.element_multiplicity(base_element_number)-first_selected_component >= n_components_, @@ -2245,24 +2214,16 @@ FEEvaluationBase active_quad_index (other.active_quad_index), matrix_info (other.matrix_info), dof_info (other.dof_info), - mapping_info (other.mapping_info), + mapping_data (other.mapping_data), data (other.matrix_info == nullptr ? new internal::MatrixFreeFunctions::ShapeInfo>(*other.data) : other.data), - cartesian_data (nullptr), jacobian (nullptr), J_value (nullptr), - quadrature_weights (mapping_info != nullptr ? - mapping_info->mapping_data_gen[quad_no]. - quadrature_weights[active_quad_index].begin() - : - nullptr), - quadrature_points (nullptr), - jacobian_grad (nullptr), - jacobian_grad_upper(nullptr), + quadrature_weights (other.matrix_info == nullptr ? nullptr : + mapping_data->descriptor[active_quad_index].quadrature_weights.begin()), cell (numbers::invalid_unsigned_int), cell_type (internal::MatrixFreeFunctions::general), - cell_data_number (numbers::invalid_unsigned_int), dof_values_initialized (false), values_quad_initialized (false), gradients_quad_initialized(false), @@ -2281,10 +2242,11 @@ FEEvaluationBase MappingDataOnTheFly(other.mapped_geometry->get_fe_values().get_mapping(), other.mapped_geometry->get_quadrature(), other.mapped_geometry->get_fe_values().get_update_flags())); - jacobian = mapped_geometry->get_inverse_jacobians().begin(); - J_value = mapped_geometry->get_JxW_values().begin(); - quadrature_points = mapped_geometry->get_quadrature_points().begin(); + mapping_data = &mapped_geometry->get_data_storage(); cell = 0; + + jacobian = mapped_geometry->get_data_storage().jacobians[0].begin(); + J_value = mapped_geometry->get_data_storage().JxW_values.begin(); } } @@ -2315,7 +2277,7 @@ FEEvaluationBase matrix_info = other.matrix_info; dof_info = other.dof_info; - mapping_info = other.mapping_info; + mapping_data = other.mapping_data; if (other.matrix_info == nullptr) { data = new internal::MatrixFreeFunctions::ShapeInfo>(*other.data); @@ -2328,20 +2290,14 @@ FEEvaluationBase } set_data_pointers(); - cartesian_data = nullptr; jacobian = nullptr; J_value = nullptr; - quadrature_weights = mapping_info != nullptr ? - mapping_info->mapping_data_gen[quad_no]. - quadrature_weights[active_quad_index].begin() - : - nullptr; - quadrature_points = nullptr; - jacobian_grad = nullptr; - jacobian_grad_upper = nullptr; + quadrature_weights = (mapping_data != nullptr ? + mapping_data->descriptor[active_quad_index].quadrature_weights.begin() + : + nullptr); cell = numbers::invalid_unsigned_int; cell_type = internal::MatrixFreeFunctions::general; - cell_data_number = numbers::invalid_unsigned_int; // Create deep copy of mapped geometry for use in parallel... if (other.mapped_geometry.get() != nullptr) @@ -2351,10 +2307,10 @@ FEEvaluationBase MappingDataOnTheFly(other.mapped_geometry->get_fe_values().get_mapping(), other.mapped_geometry->get_quadrature(), other.mapped_geometry->get_fe_values().get_update_flags())); - jacobian = mapped_geometry->get_inverse_jacobians().begin(); - J_value = mapped_geometry->get_JxW_values().begin(); - quadrature_points = mapped_geometry->get_quadrature_points().begin(); cell = 0; + mapping_data = &mapped_geometry->get_data_storage(); + jacobian = mapped_geometry->get_data_storage().jacobians[0].begin(); + J_value = mapped_geometry->get_data_storage().JxW_values.begin(); } return *this; @@ -2427,72 +2383,22 @@ FEEvaluationBase template inline void -FEEvaluationBase::reinit (const unsigned int cell_in) +FEEvaluationBase::reinit (const unsigned int cell_index) { Assert (mapped_geometry == nullptr, ExcMessage("FEEvaluation was initialized without a matrix-free object." " Integer indexing is not possible")); if (mapped_geometry != nullptr) return; - Assert (dof_info != nullptr, ExcNotInitialized()); - Assert (mapping_info != nullptr, ExcNotInitialized()); - AssertIndexRange (cell_in, dof_info->row_starts.size()-1); - AssertDimension (((dof_info->cell_active_fe_index.size() > 0) ? - dof_info->cell_active_fe_index[cell_in] : 0), - active_fe_index); - cell = cell_in; - cell_type = mapping_info->get_cell_type(cell); - cell_data_number = mapping_info->get_cell_data_index(cell); - - if (mapping_info->quadrature_points_initialized == true) - { - AssertIndexRange (cell_data_number, mapping_info-> - mapping_data_gen[quad_no].rowstart_q_points.size()); - const unsigned int index = mapping_info->mapping_data_gen[quad_no]. - rowstart_q_points[cell]; - AssertIndexRange (index, mapping_info->mapping_data_gen[quad_no]. - quadrature_points.size()); - quadrature_points = - &mapping_info->mapping_data_gen[quad_no].quadrature_points[index]; - } - if (cell_type == internal::MatrixFreeFunctions::cartesian) - { - cartesian_data = &mapping_info->cartesian_data[cell_data_number].first; - J_value = &mapping_info->cartesian_data[cell_data_number].second; - } - else if (cell_type == internal::MatrixFreeFunctions::affine) - { - jacobian = &mapping_info->affine_data[cell_data_number].first; - J_value = &mapping_info->affine_data[cell_data_number].second; - } - else - { - const unsigned int rowstart = mapping_info-> - mapping_data_gen[quad_no].rowstart_jacobians[cell_data_number]; - AssertIndexRange (rowstart, mapping_info-> - mapping_data_gen[quad_no].jacobians.size()); - jacobian = - &mapping_info->mapping_data_gen[quad_no].jacobians[rowstart]; - if (mapping_info->JxW_values_initialized == true) - { - AssertIndexRange (rowstart, mapping_info-> - mapping_data_gen[quad_no].JxW_values.size()); - J_value = &(mapping_info->mapping_data_gen[quad_no]. - JxW_values[rowstart]); - } - if (mapping_info->second_derivatives_initialized == true) - { - AssertIndexRange(rowstart, mapping_info-> - mapping_data_gen[quad_no].jacobians_grad_diag.size()); - jacobian_grad = &mapping_info->mapping_data_gen[quad_no]. - jacobians_grad_diag[rowstart]; - AssertIndexRange(rowstart, mapping_info-> - mapping_data_gen[quad_no].jacobians_grad_upper.size()); - jacobian_grad_upper = &mapping_info->mapping_data_gen[quad_no]. - jacobians_grad_upper[rowstart]; - } - } + Assert (this->dof_info != nullptr, ExcNotInitialized()); + Assert (this->mapping_data != nullptr, ExcNotInitialized()); + this->cell = cell_index; + this->cell_type = this->matrix_info->get_mapping_info().get_cell_type(cell_index); + + const unsigned int offsets = this->mapping_data->data_index_offsets[cell_index]; + this->jacobian = &this->mapping_data->jacobians[0][offsets]; + this->J_value = &this->mapping_data->JxW_values[offsets]; #ifdef DEBUG dof_values_initialized = false; @@ -2550,8 +2456,24 @@ unsigned int FEEvaluationBase ::get_cell_data_number () const { - Assert (cell != numbers::invalid_unsigned_int, ExcNotInitialized()); - return cell_data_number; + return get_mapping_data_index_offset(); +} + + + +template +inline +unsigned int +FEEvaluationBase +::get_mapping_data_index_offset () const +{ + if (matrix_info == 0) + return 0; + else + { + AssertIndexRange(cell, this->mapping_data->data_index_offsets.size()); + return this->mapping_data->data_index_offsets[cell]; + } } @@ -2589,7 +2511,7 @@ FEEvaluationBase if (this->cell_type == internal::MatrixFreeFunctions::cartesian || this->cell_type == internal::MatrixFreeFunctions::affine) { - Assert (this->mapping_info != nullptr, ExcNotImplemented()); + Assert (this->mapping_data != nullptr, ExcNotImplemented()); VectorizedArray J = this->J_value[0]; for (unsigned int q=0; qdata->n_q_points; ++q) JxW_values[q] = J * this->quadrature_weights[q]; @@ -2606,11 +2528,12 @@ inline VectorizedArray FEEvaluationBase::JxW(const unsigned int q_point) const { + AssertIndexRange(q_point, data->n_q_points); Assert (this->J_value != nullptr, ExcNotInitialized()); if (this->cell_type == internal::MatrixFreeFunctions::cartesian || this->cell_type == internal::MatrixFreeFunctions::affine) { - Assert (this->mapping_info != nullptr, ExcInternalError()); + Assert (this->quadrature_weights != nullptr, ExcInternalError()); return this->J_value[0] * this->quadrature_weights[q_point]; } else @@ -2619,6 +2542,22 @@ FEEvaluationBase::JxW(const unsigned int q_point) cons +template +inline +Tensor<2,dim,VectorizedArray > +FEEvaluationBase +::inverse_jacobian(const unsigned int q_index) const +{ + AssertIndexRange(q_index, data->n_q_points); + Assert (this->jacobian != nullptr, ExcNotImplemented()); + if (this->cell_type <= internal::MatrixFreeFunctions::affine) + return jacobian[0]; + else + return jacobian[q_index]; +} + + + namespace internal { // write access to generic vectors that have operator (). @@ -3749,13 +3688,13 @@ FEEvaluationBase for (unsigned int comp=0; compgradients_quad[comp][d][q_point] * - cartesian_data[0][d]); + jacobian[0][d][d]); } // cell with general/affine Jacobian else { const Tensor<2,dim,VectorizedArray > &jac = - this->cell_type == internal::MatrixFreeFunctions::general ? + this->cell_type > internal::MatrixFreeFunctions::affine ? jacobian[q_point] : jacobian[0]; for (unsigned int comp=0; comp internal::ExcAccessToUninitializedField()); AssertIndexRange (q_point, this->data->n_q_points); + Assert(jacobian != nullptr, ExcNotImplemented()); + const Tensor<2,dim,VectorizedArray > &jac = + jacobian[this->cell_type <= internal::MatrixFreeFunctions::affine ? 0 : q_point]; + Tensor<2,dim,VectorizedArray > hessian_out [n_components]; // Cartesian cell if (this->cell_type == internal::MatrixFreeFunctions::cartesian) { - const Tensor<1,dim,VectorizedArray > &jac = cartesian_data[0]; for (unsigned int comp=0; comphessians_quad[comp][d][q_point] * - jac[d] * jac[d]); + jac[d][d] * jac[d][d]); switch (dim) { case 1: break; case 2: hessian_out[comp][0][1] = (this->hessians_quad[comp][2][q_point] * - jac[0] * jac[1]); + jac[0][0] * jac[1][1]); break; case 3: hessian_out[comp][0][1] = (this->hessians_quad[comp][3][q_point] * - jac[0] * jac[1]); + jac[0][0] * jac[1][1]); hessian_out[comp][0][2] = (this->hessians_quad[comp][4][q_point] * - jac[0] * jac[2]); + jac[0][0] * jac[2][2]); hessian_out[comp][1][2] = (this->hessians_quad[comp][5][q_point] * - jac[1] * jac[2]); + jac[1][1] * jac[2][2]); break; default: Assert (false, ExcNotImplemented()); @@ -3875,16 +3817,9 @@ FEEvaluationBase hessian_out[comp][e][d] = hessian_out[comp][d][e]; } } - // cell with general Jacobian - else if (this->cell_type == internal::MatrixFreeFunctions::general) + // cell with general Jacobian, but constant within the cell + else if (this->cell_type == internal::MatrixFreeFunctions::affine) { - Assert (this->mapping_info->second_derivatives_initialized == true, - ExcNotInitialized()); - const Tensor<2,dim,VectorizedArray > &jac = jacobian[q_point]; - const Tensor<2,dim,VectorizedArray > &jac_grad = jacobian_grad[q_point]; - const Tensor<1,(dim>1?dim*(dim-1)/2:1), - Tensor<1,dim,VectorizedArray > > - & jac_grad_UT = jacobian_grad_upper[q_point]; for (unsigned int comp=0; comp hessian_out[comp][d][e] += jac[d][f] * tmp[f][e]; } - // add diagonal part of J' * grad(u) - for (unsigned int d=0; dgradients_quad[comp][e][q_point]); - - // add off-diagonal part of J' * grad(u) - for (unsigned int d=0, count=0; dgradients_quad[comp][f][q_point]); + // no J' * grad(u) part here because the Jacobian is constant + // throughout the cell and hence, its derivative is zero // take symmetric part for (unsigned int d=0; d hessian_out[comp][e][d] = hessian_out[comp][d][e]; } } - // cell with general Jacobian, but constant within the cell - else // if (this->cell_type == internal::MatrixFreeFunctions::affine) + // cell with general Jacobian + else { - const Tensor<2,dim,VectorizedArray > &jac = jacobian[0]; + const Tensor<1,dim*(dim+1)/2,Tensor<1,dim,VectorizedArray > > &jac_grad = + mapping_data->jacobian_gradients[0][this->get_mapping_data_index_offset()+q_point]; for (unsigned int comp=0; comp hessian_out[comp][d][e] += jac[d][f] * tmp[f][e]; } - // no J' * grad(u) part here because the Jacobian is constant - // throughout the cell and hence, its derivative is zero + // add diagonal part of J' * grad(u) + for (unsigned int d=0; dgradients_quad[comp][e][q_point]); + + // add off-diagonal part of J' * grad(u) + for (unsigned int d=0, count=dim; dgradients_quad[comp][f][q_point]); // take symmetric part for (unsigned int d=0; d internal::ExcAccessToUninitializedField()); AssertIndexRange (q_point, this->data->n_q_points); + Assert(jacobian != nullptr, ExcNotImplemented()); + const Tensor<2,dim,VectorizedArray > &jac = + jacobian[this->cell_type <= internal::MatrixFreeFunctions::affine ? 0 : q_point]; + Tensor<1,n_components_,Tensor<1,dim,VectorizedArray > > hessian_out; // Cartesian cell if (this->cell_type == internal::MatrixFreeFunctions::cartesian) { - const Tensor<1,dim,VectorizedArray > &jac = cartesian_data[0]; for (unsigned int comp=0; comphessians_quad[comp][d][q_point] * - jac[d] * jac[d]); + jac[d][d] * jac[d][d]); } - // cell with general Jacobian - else if (this->cell_type == internal::MatrixFreeFunctions::general) + // cell with general Jacobian, but constant within the cell + else if (this->cell_type == internal::MatrixFreeFunctions::affine) { - Assert (this->mapping_info->second_derivatives_initialized == true, - ExcNotInitialized()); - const Tensor<2,dim,VectorizedArray > &jac = jacobian[q_point]; - const Tensor<2,dim,VectorizedArray > &jac_grad = jacobian_grad[q_point]; for (unsigned int comp=0; comp for (unsigned int f=1; fgradients_quad[comp][e][q_point]); } } - // cell with general Jacobian, but constant within the cell - else // if (this->cell_type == internal::MatrixFreeFunctions::affine) + // cell with general Jacobian + else { - const Tensor<2,dim,VectorizedArray > &jac = jacobian[0]; + const Tensor<1,dim*(dim+1)/2,Tensor<1,dim,VectorizedArray > > &jac_grad = + mapping_data->jacobian_gradients[0][this->get_mapping_data_index_offset()+q_point]; for (unsigned int comp=0; comp for (unsigned int f=1; fgradients_quad[comp][e][q_point]); } } return hessian_out; @@ -4114,22 +4050,25 @@ FEEvaluationBase Assert (this->cell != numbers::invalid_unsigned_int, ExcNotInitialized()); AssertIndexRange (q_point, this->data->n_q_points); this->gradients_quad_submitted = true; + Assert (this->J_value != nullptr, ExcNotInitialized()); + Assert (this->jacobian != nullptr, ExcNotInitialized()); #endif + if (this->cell_type == internal::MatrixFreeFunctions::cartesian) { const VectorizedArray JxW = J_value[0] * quadrature_weights[q_point]; for (unsigned int comp=0; compgradients_quad[comp][d][q_point] = (grad_in[comp][d] * - cartesian_data[0][d] * JxW); + jacobian[0][d][d] * JxW); } else { const Tensor<2,dim,VectorizedArray > &jac = - this->cell_type == internal::MatrixFreeFunctions::general ? + this->cell_type > internal::MatrixFreeFunctions::affine ? jacobian[q_point] : jacobian[0]; const VectorizedArray JxW = - this->cell_type == internal::MatrixFreeFunctions::general ? + this->cell_type > internal::MatrixFreeFunctions::affine ? J_value[q_point] : J_value[0] * quadrature_weights[q_point]; for (unsigned int comp=0; comp internal::ExcAccessToUninitializedField()); AssertIndexRange (q_point, this->data->n_q_points); + Assert (this->jacobian != nullptr, ExcNotInitialized()); + Tensor<1,dim,VectorizedArray > grad_out; - // Cartesian cell if (this->cell_type == internal::MatrixFreeFunctions::cartesian) { for (unsigned int d=0; dgradients_quad[0][d][q_point] * - this->cartesian_data[0][d]); + this->jacobian[0][d][d]); } - // cell with general/constant Jacobian + // cell with general/affine Jacobian else { const Tensor<2,dim,VectorizedArray > &jac = - this->cell_type == internal::MatrixFreeFunctions::general ? - this->jacobian[q_point] : this->jacobian[0]; + this->jacobian[this->cell_type > internal::MatrixFreeFunctions::affine ? + q_point : 0]; for (unsigned int d=0; dgradients_quad[0][0][q_point]); + grad_out[d] = jac[d][0] * this->gradients_quad[0][0][q_point]; for (unsigned int e=1; egradients_quad[0][e][q_point]); + grad_out[d] += jac[d][e] * this->gradients_quad[0][e][q_point]; } } return grad_out; @@ -4427,36 +4367,40 @@ inline void FEEvaluationAccess ::submit_gradient (const Tensor<1,dim,VectorizedArray > grad_in, - const unsigned int q_point) + const unsigned int q_index) { #ifdef DEBUG Assert (this->cell != numbers::invalid_unsigned_int, ExcNotInitialized()); - AssertIndexRange (q_point, this->data->n_q_points); + AssertIndexRange (q_index, this->data->n_q_points); this->gradients_quad_submitted = true; + Assert (this->J_value != nullptr, ExcNotInitialized()); + Assert (this->jacobian != nullptr, ExcNotInitialized()); #endif + if (this->cell_type == internal::MatrixFreeFunctions::cartesian) { - const VectorizedArray JxW = this->J_value[0] * this->quadrature_weights[q_point]; + const VectorizedArray JxW = this->J_value[0] * + this->quadrature_weights[q_index]; for (unsigned int d=0; dgradients_quad[0][d][q_point] = (grad_in[d] * - this->cartesian_data[0][d] * + this->gradients_quad[0][d][q_index] = (grad_in[d] * + this->jacobian[0][d][d] * JxW); } // general/affine cell type else { const Tensor<2,dim,VectorizedArray > &jac = - this->cell_type == internal::MatrixFreeFunctions::general ? - this->jacobian[q_point] : this->jacobian[0]; + this->cell_type > internal::MatrixFreeFunctions::affine ? + this->jacobian[q_index] : this->jacobian[0]; const VectorizedArray JxW = - this->cell_type == internal::MatrixFreeFunctions::general ? - this->J_value[q_point] : this->J_value[0] * this->quadrature_weights[q_point]; + this->cell_type > internal::MatrixFreeFunctions::affine ? + this->J_value[q_index] : this->J_value[0] * this->quadrature_weights[q_index]; for (unsigned int d=0; d new_val = jac[0][d] * grad_in[0]; for (unsigned int e=1; egradients_quad[0][d][q_point] = new_val * JxW; + this->gradients_quad[0][d][q_index] = new_val * JxW; } } } @@ -4552,6 +4496,7 @@ FEEvaluationAccess Assert (this->gradients_quad_initialized==true, internal::ExcAccessToUninitializedField()); AssertIndexRange (q_point, this->data->n_q_points); + Assert (this->jacobian != nullptr, ExcNotInitialized()); VectorizedArray divergence; @@ -4559,10 +4504,10 @@ FEEvaluationAccess if (this->cell_type == internal::MatrixFreeFunctions::cartesian) { divergence = (this->gradients_quad[0][0][q_point] * - this->cartesian_data[0][0]); + this->jacobian[0][0][0]); for (unsigned int d=1; dgradients_quad[d][d][q_point] * - this->cartesian_data[0][d]); + this->jacobian[0][d][d]); } // cell with general/constant Jacobian else @@ -4714,7 +4659,10 @@ FEEvaluationAccess Assert (this->cell != numbers::invalid_unsigned_int, ExcNotInitialized()); AssertIndexRange (q_point, this->data->n_q_points); this->gradients_quad_submitted = true; + Assert (this->J_value != nullptr, ExcNotInitialized()); + Assert (this->jacobian != nullptr, ExcNotInitialized()); #endif + if (this->cell_type == internal::MatrixFreeFunctions::cartesian) { const VectorizedArray fac = this->J_value[0] * @@ -4722,7 +4670,7 @@ FEEvaluationAccess for (unsigned int d=0; dgradients_quad[d][d][q_point] = (fac * - this->cartesian_data[0][d]); + this->jacobian[0][d][d]); for (unsigned int e=d+1; egradients_quad[d][e][q_point] = VectorizedArray(); @@ -4764,22 +4712,25 @@ FEEvaluationAccess Assert (this->cell != numbers::invalid_unsigned_int, ExcNotInitialized()); AssertIndexRange (q_point, this->data->n_q_points); this->gradients_quad_submitted = true; + Assert (this->J_value != nullptr, ExcNotInitialized()); + Assert (this->jacobian != nullptr, ExcNotInitialized()); #endif + if (this->cell_type == internal::MatrixFreeFunctions::cartesian) { const VectorizedArray JxW = this->J_value[0] * this->quadrature_weights[q_point]; for (unsigned int d=0; dgradients_quad[d][d][q_point] = (sym_grad.access_raw_entry(d) * JxW * - this->cartesian_data[0][d]); + this->jacobian[0][d][d]); for (unsigned int e=0, counter=dim; e value = sym_grad.access_raw_entry(counter) * JxW; this->gradients_quad[e][d][q_point] = (value * - this->cartesian_data[0][d]); + this->jacobian[0][d][d]); this->gradients_quad[d][e][q_point] = (value * - this->cartesian_data[0][e]); + this->jacobian[0][e][e]); } } // general/affine cell type @@ -4943,23 +4894,13 @@ FEEvaluationAccess<1,1,Number> internal::ExcAccessToUninitializedField()); AssertIndexRange (q_point, this->data->n_q_points); - Tensor<1,1,VectorizedArray > grad_out; + const Tensor<2,1,VectorizedArray > &jac = + this->cell_type == internal::MatrixFreeFunctions::general ? + this->jacobian[q_point] : this->jacobian[0]; - // Cartesian cell - if (this->cell_type == internal::MatrixFreeFunctions::cartesian) - { - grad_out[0] = (this->gradients_quad[0][0][q_point] * - this->cartesian_data[0][0]); - } - // cell with general/constant Jacobian - else - { - const Tensor<2,1,VectorizedArray > &jac = - this->cell_type == internal::MatrixFreeFunctions::general ? - this->jacobian[q_point] : this->jacobian[0]; + Tensor<1,1,VectorizedArray > grad_out; + grad_out[0] = jac[0][0] * this->gradients_quad[0][0][q_point]; - grad_out[0] = (jac[0][0] * this->gradients_quad[0][0][q_point]); - } return grad_out; } @@ -5052,25 +4993,15 @@ FEEvaluationAccess<1,1,Number> AssertIndexRange (q_point, this->data->n_q_points); this->gradients_quad_submitted = true; #endif - if (this->cell_type == internal::MatrixFreeFunctions::cartesian) - { - const VectorizedArray JxW = this->J_value[0] * this->quadrature_weights[q_point]; - this->gradients_quad[0][0][q_point] = (grad_in[0] * - this->cartesian_data[0][0] * - JxW); - } - // general/affine cell type - else - { - const Tensor<2,1,VectorizedArray > &jac = - this->cell_type == internal::MatrixFreeFunctions::general ? - this->jacobian[q_point] : this->jacobian[0]; - const VectorizedArray JxW = - this->cell_type == internal::MatrixFreeFunctions::general ? - this->J_value[q_point] : this->J_value[0] * this->quadrature_weights[q_point]; - this->gradients_quad[0][0][q_point] = jac[0][0] * grad_in[0] * JxW; - } + const Tensor<2,1,VectorizedArray > &jac = + this->cell_type > internal::MatrixFreeFunctions::affine ? + this->jacobian[q_point] : this->jacobian[0]; + const VectorizedArray JxW = + this->cell_type > internal::MatrixFreeFunctions::affine ? + this->J_value[q_point] : this->J_value[0] * this->quadrature_weights[q_point]; + + this->gradients_quad[0][0][q_point] = jac[0][0] * grad_in[0] * JxW; } @@ -5252,13 +5183,13 @@ FEEvaluation proposed_dof_comp = no; break; } - if (static_n_q_points == - this->mapping_info->mapping_data_gen[this->quad_no].n_q_points[this->active_quad_index]) + if (n_q_points == + this->mapping_data->descriptor[this->active_quad_index].n_q_points) proposed_quad_comp = this->quad_no; else - for (unsigned int no=0; nomapping_info->mapping_data_gen.size(); ++no) - if (this->mapping_info->mapping_data_gen[no].n_q_points[this->active_quad_index] - == static_n_q_points) + for (unsigned int no=0; nomatrix_info->get_mapping_info().cell_data.size(); ++no) + if (this->matrix_info->get_mapping_info().cell_data[no].descriptor[this->active_quad_index].n_q_points + == n_q_points) { proposed_quad_comp = no; break; @@ -5326,8 +5257,7 @@ FEEvaluation if (fe_no != numbers::invalid_unsigned_int) { AssertDimension (n_q_points, - this->mapping_info->mapping_data_gen[this->quad_no]. - n_q_points[this->active_quad_index]); + this->mapping_data->descriptor[this->active_quad_index].n_q_points); AssertDimension (this->data->dofs_per_component_on_cell * this->n_fe_components, this->dof_info->dofs_per_cell[this->active_fe_index]); } @@ -5343,10 +5273,24 @@ Point > FEEvaluation ::quadrature_point (const unsigned int q) const { - Assert (this->mapping_info->quadrature_points_initialized == true, - ExcNotInitialized()); - Assert (this->quadrature_points != nullptr, ExcNotInitialized()); + if (this->matrix_info == nullptr) + { + Assert ((this->mapped_geometry->get_fe_values().get_update_flags() | + update_quadrature_points), + ExcNotInitialized()); + } + else + { + Assert(this->mapping_data->quadrature_point_offsets.empty() == false, + ExcNotInitialized()); + } + AssertIndexRange (q, n_q_points); + const Point > *quadrature_points = &this->mapping_data-> + quadrature_points[this->mapping_data->quadrature_point_offsets[this->cell]]; + + const unsigned int n_q_points_1d_actual = + fe_degree == -1 ? this->data->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 @@ -5357,15 +5301,15 @@ FEEvaluation switch (dim) { case 1: - return this->quadrature_points[q]; + return quadrature_points[q]; case 2: - point[0] = this->quadrature_points[q%n_q_points_1d][0]; - point[1] = this->quadrature_points[q/n_q_points_1d][1]; + 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] = this->quadrature_points[q%n_q_points_1d][0]; - point[1] = this->quadrature_points[(q/n_q_points_1d)%n_q_points_1d][1]; - point[2] = this->quadrature_points[q/(n_q_points_1d*n_q_points_1d)][2]; + 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()); @@ -5374,7 +5318,7 @@ FEEvaluation } // all other cases: just return the respective data as it is fully stored else - return this->quadrature_points[q]; + return quadrature_points[q]; } diff --git a/include/deal.II/matrix_free/helper_functions.h b/include/deal.II/matrix_free/helper_functions.h index 8cb01efebb..1720820106 100644 --- a/include/deal.II/matrix_free/helper_functions.h +++ b/include/deal.II/matrix_free/helper_functions.h @@ -132,64 +132,6 @@ namespace internal unsigned int n_procs; }; - /** - * Data type to identify cell type. - */ - enum CellType - { - /** - * The cell is Cartesian. - */ - cartesian = 0, - /** - * The cell may be described with an affine mapping. - */ - affine = 1, - /** - * There is no special information available for compressing the - * representation of the cell. - */ - general = 2, - /** - * The cell type is undefined. - */ - undefined = 3 - }; - - - /** - * A class that is used to compare floating point arrays (e.g. - * std::vectors, Tensor<1,dim>, etc.). The idea of this class is to - * consider two arrays as equal if they are the same within a given - * tolerance. We use this comparator class within an 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 - struct FPArrayComparator - { - FPArrayComparator (const Number scaling); - - bool operator() (const std::vector &v1, - const std::vector &v2) const; - - template - bool operator ()(const Tensor<1,dim,Tensor<1,VectorizedArray::n_array_elements,Number> > &t1, - const Tensor<1,dim,Tensor<1,VectorizedArray::n_array_elements,Number> > &t2) const; - - template - bool operator ()(const Tensor<2,dim,Tensor<1,VectorizedArray::n_array_elements,Number> > &t1, - const Tensor<2,dim,Tensor<1,VectorizedArray::n_array_elements,Number> > &t2) const; - - Number tolerance; - }; - - // Note: Implementation in matrix_free.templates.h - } // end of namespace MatrixFreeFunctions } // end of namespace internal diff --git a/include/deal.II/matrix_free/mapping_data_on_the_fly.h b/include/deal.II/matrix_free/mapping_data_on_the_fly.h index a5c196a605..265731fae6 100644 --- a/include/deal.II/matrix_free/mapping_data_on_the_fly.h +++ b/include/deal.II/matrix_free/mapping_data_on_the_fly.h @@ -27,6 +27,7 @@ #include #include #include +#include DEAL_II_NAMESPACE_OPEN @@ -101,36 +102,13 @@ namespace internal const dealii::FEValues &get_fe_values () const; /** - * Return a vector of inverse transpose Jacobians. For compatibility - * with FEEvaluation, it returns tensors of vectorized arrays, even - * though all components are equal. + * Return a reference to the underlying storage field of type + * MappingInfoStorage of the same format as the data fields in + * MappingInfo. This ensures compatibility with the precomputed data + * fields in the MappingInfo class. */ - const AlignedVector > > & - get_inverse_jacobians() const; - - /** - * Return a vector of quadrature weights times the Jacobian determinant - * (JxW). For compatibility with FEEvaluation, it returns tensors of - * vectorized arrays, even though all components are equal. - */ - const AlignedVector > & - get_JxW_values() const; - - /** - * Return a vector of quadrature points in real space on the given cell. - * For compatibility with FEEvaluation, it returns tensors of vectorized - * arrays, even though all components are equal. - */ - const AlignedVector > > & - get_quadrature_points() const; - - /** - * Return a vector of normal vectors in real space on the given cell. - * For compatibility with FEEvaluation, it returns tensors of vectorized - * arrays, even though all components are equal. - */ - const AlignedVector > > & - get_normal_vectors() const; + const MappingInfoStorage & + get_data_storage() const; /** * Return a reference to 1D quadrature underlying this object. @@ -162,24 +140,10 @@ namespace internal const Quadrature<1> quadrature_1d; /** - * Inverse Jacobians, stored in vectorized array form. - */ - AlignedVector > > inverse_jacobians; - - /** - * Stored Jacobian determinants and quadrature weights - */ - AlignedVector > jxw_values; - - /** - * Stored quadrature points + * The storage part created for a single cell and held in analogy to + * MappingInfo. */ - AlignedVector > > quadrature_points; - - /** - * Stored normal vectors (for face integration) - */ - AlignedVector > > normal_vectors; + MappingInfoStorage mapping_info_storage; }; @@ -192,13 +156,24 @@ namespace internal const UpdateFlags update_flags) : fe_values(mapping, fe_dummy, Quadrature(quadrature), - internal::MatrixFreeFunctions::MappingInfo::compute_update_flags(update_flags)), - quadrature_1d(quadrature), - inverse_jacobians(fe_values.get_quadrature().size()), - jxw_values(fe_values.get_quadrature().size()), - quadrature_points(fe_values.get_quadrature().size()), - normal_vectors(fe_values.get_quadrature().size()) + MappingInfo::compute_update_flags(update_flags)), + quadrature_1d(quadrature) { + mapping_info_storage.descriptor.resize(1); + mapping_info_storage.descriptor[0].initialize(quadrature); + mapping_info_storage.data_index_offsets.resize(1); + mapping_info_storage.JxW_values.resize(fe_values.n_quadrature_points); + mapping_info_storage.jacobians[0].resize(fe_values.n_quadrature_points); + if (update_flags & update_quadrature_points) + { + mapping_info_storage.quadrature_point_offsets.resize(1, 0); + mapping_info_storage.quadrature_points.resize(fe_values.n_quadrature_points); + } + if (fe_values.get_update_flags() & update_normal_vectors) + { + mapping_info_storage.normal_vectors.resize(fe_values.n_quadrature_points); + mapping_info_storage.normals_times_jacobians[0].resize(fe_values.n_quadrature_points); + } Assert(!(fe_values.get_update_flags() & update_jacobian_grads), ExcNotImplemented()); } @@ -210,17 +185,9 @@ namespace internal MappingDataOnTheFly::MappingDataOnTheFly (const Quadrature<1> &quadrature, const UpdateFlags update_flags) : - fe_values(fe_dummy, Quadrature(quadrature), - internal::MatrixFreeFunctions::MappingInfo::compute_update_flags(update_flags)), - quadrature_1d(quadrature), - inverse_jacobians(fe_values.get_quadrature().size()), - jxw_values(fe_values.get_quadrature().size()), - quadrature_points(fe_values.get_quadrature().size()), - normal_vectors(fe_values.get_quadrature().size()) - { - Assert(!(fe_values.get_update_flags() & update_jacobian_grads), - ExcNotImplemented()); - } + MappingDataOnTheFly(::dealii::StaticMappingQ1::mapping, + quadrature, update_flags) + {} @@ -235,18 +202,27 @@ namespace internal fe_values.reinit(present_cell); for (unsigned int q=0; q jac = fe_values.jacobian(q); + jac = invert(transpose(jac)); + for (unsigned int d=0; d inline - const AlignedVector > > & - MappingDataOnTheFly::get_inverse_jacobians() const - { - return inverse_jacobians; - } - - - - template - inline - const AlignedVector > > & - MappingDataOnTheFly::get_normal_vectors() const - { - return normal_vectors; - } - - - - template - inline - const AlignedVector > > & - MappingDataOnTheFly::get_quadrature_points() const - { - return quadrature_points; - } - - - - template - inline - const AlignedVector > & - MappingDataOnTheFly::get_JxW_values() const + const MappingInfoStorage & + MappingDataOnTheFly::get_data_storage() const { - return jxw_values; + return mapping_info_storage; } diff --git a/include/deal.II/matrix_free/mapping_info.h b/include/deal.II/matrix_free/mapping_info.h index d4afeaea8a..74da08ff99 100644 --- a/include/deal.II/matrix_free/mapping_info.h +++ b/include/deal.II/matrix_free/mapping_info.h @@ -37,71 +37,268 @@ namespace internal namespace MatrixFreeFunctions { /** - * The class that stores all geometry-dependent data related with cell - * interiors for use in the matrix-free class. + * An enum to identify various types of cells. The most general type is + * what we typically compute in the FEValues context but for many + * geometries we can save significant storage. + */ + enum CellType : unsigned char + { + /** + * The cell is Cartesian. + */ + cartesian = 0, + /** + * The cell may be described with an affine mapping. + */ + affine = 1, + /** + * The current face is flat, i.e., the normal factor on a face is the + * same on all quadrature points. + */ + flat_faces = 2, + /** + * There is no special information available for compressing the + * representation of the cell. + */ + general = 3 + }; + + + + /** + * Definition of a structure that stores all cached data related to the + * evaluated geometry from the mapping. In order to support hp-adaptivity + * and compressed storage (in particular for Jacobians, JxW values, and + * normals), storage length can be different for different rows. Thus, it + * allows to jump at the data of individual rows similar to compressed row + * storage in sparse matrices. We have two different start indices for + * fields with different sizes. The first category of offsets are the + * indices for Jacobians of the transformation from unit to real cell (we + * store the inverse Jacobian), second derivatives, JxW values, and normal + * vectors. We keep separate arrays for all these data structures because + * a user code might access only some of them. In such a case, one array + * will be gone through in a contiguous order with access to all entries, + * which makes it easy for the processor to prefetch data. Having all data + * in a single array would require some strides in the access pattern, + * which is much more complicated for the processor to predict (and indeed + * leads to prefetching of data that does not get used on Intel processors + * such as BroadwellEP). * - * @author Katharina Kormann and Martin Kronbichler, 2010, 2011 + * The second category of indices are the offsets for the quadrature + * points. Quadrature points can be compressed less than the other fields + * and thus need longer fields. Quadrature point indices are often used in + * other contexts such as evaluation of right hand sides. + * + * The third component is a descriptor of data from the unit cells, called + * QuadratureDescriptor, which contains the quadrature weights and + * permutations of how to go through quadrature points in case of face + * data. The latter comes in a vector for the support of hp adaptivity, + * with several data fields for the individual quadrature formulas. + * + * @author Katharina Kormann, Martin Kronbichler, 2018 */ - template - struct MappingInfo + template + struct MappingInfoStorage { + struct QuadratureDescriptor + { + /** + * Constructor. Does nothing. + */ + QuadratureDescriptor(); + + /** + * Set up the lengths in the various members of this struct. + */ + void initialize(const Quadrature<1> &quadrature_1d, + const UpdateFlags update_flags_inner_faces = update_default); + + /** + * Returns the memory consumption in bytes. + */ + std::size_t memory_consumption() const; + + /** + * Number of quadrature points applied on the given cell or face. + */ + unsigned int n_q_points; + + /** + * Quadrature formula applied on the given cell or face. + */ + Quadrature quadrature; + + /** + * Quadrature weights separated by dimension for use in specific + * situations. + */ + std::array, structdim> tensor_quadrature_weights; + + /** + * A cached vector of quadrature weights in the given number format + * (non-vectorized, as it is cheap to broadcast the value to all lanes + * when it is used in a vectorized context). + */ + AlignedVector quadrature_weights; + + /** + * For quadrature on faces, the evaluation of basis functions is not + * in the correct order if a face is not in the standard orientation + * to a given element. This data structure is used to re-order the + * data evaluated on quadrature points to represent the correct order. + */ + dealii::Table<2,unsigned int> face_orientations; + }; + /** - * Determines how many bits of an unsigned int are used to distinguish - * the cell types (Cartesian, with constant Jacobian, or general) + * A class describing the layout of the sections in the @p data_storage + * field and also includes some data that depends on the number of + * quadrature points in the hp context such as the inner quadrature + * formula and re-indexing for faces that are not in the standard + * orientation. */ - static const std::size_t n_cell_type_bits = 2; + std::vector descriptor; /** - * Determines how many types of different cells can be detected at most. - * Corresponds to the number of bits we reserved for it. + * Stores the index offset into the arrays @p jxw_values, @p jacobians, + * @p normal_vectors and the second derivatives. Note that affine cells + * have shorter fields of length 1, where the others have lengths equal + * to the numer of quadrature points of the given cell. */ - static const unsigned int n_cell_types = 1U< data_index_offsets; /** - * An abbreviation for the length of vector lines of the current data - * type. + * The storage of the Jacobian determinant (times the quadrature weight + * in case the transformation is non-affine) on quadrature + * points. + * + * Indexed by @p data_index_offsets. */ - static const unsigned int n_vector_elements = VectorizedArray::n_array_elements; + AlignedVector > JxW_values; + /** + * Stores the normal vectors. + * + * Indexed by @p data_index_offsets. + */ + AlignedVector > > normal_vectors; + + /** + * The storage of covariant transformation on quadrature points, i.e., + * the inverse and transposed Jacobians of the transformation from the + * unit to the real cell. + * + * Indexed by @p data_index_offsets. + * + * Contains two fields for access from both sides for interior faces, + * but the default case (cell integrals or boundary integrals) only + * fills the zeroth component and ignores the first one. + */ + AlignedVector > > jacobians[2]; + + /** + * The storage of the gradients of the inverse Jacobian + * transformation. Because of symmetry, only the upper diagonal and + * diagonal part are needed. The first index runs through the + * derivatives, starting with the diagonal and then continuing row-wise, + * i.e., $\partial^2/\partial x_1 \partial x_2$ first, then + * $\partial^2/\partial x_1 \partial x_3$, and so on. The second index + * is the spatial coordinate. + * + * Indexed by @p data_index_offsets. + * + * Contains two fields for access from both sides for interior faces, + * but the default case (cell integrals or boundary integrals) only + * fills the zeroth component and ignores the first one. + */ + AlignedVector > > > jacobian_gradients[2]; + + /** + * Stores the Jacobian transformations times the normal vector (this + * represents a shortcut that is accessed often and can thus get higher + * performance). + * + * Indexed by @p data_index_offsets. + */ + AlignedVector > > normals_times_jacobians [2]; + + /** + * Stores the index offset of a particular cell into the quadrature + * points array in real coordinates. Note that Cartesian cells have + * shorter fields (length is @p n_q_points_1d) than non-Cartesian cells + * (length is @p n_q_points) or faces. + */ + AlignedVector quadrature_point_offsets; + + /** + * Stores the quadrature points in real coordinates, including a + * compression scheme for Cartesian cells where we do not need to store + * the full data on all points. + * + * Indexed by @p quadrature_point_offsets. + */ + AlignedVector > > quadrature_points; + + /** + * Returns the quadrature index for a given number of quadrature + * points. If not in hp mode or if the index is not found, this + * function always returns index 0. Hence, this function does not + * check whether the given degree is actually present. + */ + unsigned int + quad_index_from_n_q_points (const unsigned int n_q_points) const; + + /** + * Prints a detailed summary of memory consumption in the different + * structures of this class to the given output stream. + */ + template + void print_memory_consumption(StreamType &out, + const SizeInfo &task_info) const; + + /** + * Returns the memory consumption in bytes. + */ + std::size_t memory_consumption () const; + }; + + + + /** + * The class that stores all geometry-dependent data related with cell + * interiors for use in the matrix-free class. + * + * @author Katharina Kormann and Martin Kronbichler, 2010, 2011, 2017 + */ + template + struct MappingInfo + { /** * Empty constructor. */ MappingInfo(); /** - * Compute the information in the given cells. The cells are specified - * by the level and the index within the level (as given by + * Compute the information in the given cells and faces. The cells are + * specified by the level and the index within the level (as given by * CellIterator::level() and CellIterator::index(), in order to allow * for different kinds of iterators, e.g. standard DoFHandler, * multigrid, etc.) on a fixed Triangulation. In addition, a mapping * and several quadrature formulas are given. */ void initialize (const dealii::Triangulation &tria, - const std::vector > &cells, + const std::vector > &cells, const std::vector &active_fe_index, const Mapping &mapping, const std::vector > &quad, - const UpdateFlags update_flags); - - /** - * Helper function to determine which update flags must be set in the - * internal functions to initialize all data as requested by the user. - */ - static UpdateFlags - compute_update_flags (const UpdateFlags update_flags, - const std::vector > &quad = - std::vector >()); + const UpdateFlags update_flags_cells); /** * Return the type of a given cell as detected during initialization. */ CellType get_cell_type (const unsigned int cell_chunk_no) const; - /** - * Return the type of a given cell as detected during initialization. - */ - unsigned int get_cell_data_index (const unsigned int cell_chunk_no) const; - /** * Clear all data fields in this class. */ @@ -118,232 +315,143 @@ namespace internal */ template void print_memory_consumption(StreamType &out, - const SizeInfo &size_info) const; + const SizeInfo &task_info) const; /** - * Stores whether a cell is Cartesian, has constant transform data - * (Jacobians) or is general. cell_type % 4 gives this information (0: - * Cartesian, 1: constant Jacobian throughout cell, 2: general cell), - * and cell_type / 4 gives the index in the data field of where to find - * the information in the fields Jacobian and JxW values (except for - * quadrature points, for which the index runs as usual). + * Stores whether a cell is Cartesian (cell type 0), has constant + * transform data (Jacobians) (cell type 1), or is general (cell type + * 3). Type 2 is only used for faces and no cells are assigned this + * value. */ - std::vector cell_type; + std::vector cell_type; /** - * The first field stores the inverse Jacobian for Cartesian cells: - * There, it is a diagonal rank-2 tensor, so we actually just store a - * rank-1 tensor. It is the same on all cells, therefore we only store - * it once per cell, and use similarities from one cell to another, too - * (on structured meshes, there are usually many cells with the same - * Jacobian). - * - * The second field stores the Jacobian determinant for Cartesian cells - * (without the quadrature weight, which depends on the quadrature - * point, whereas the determinant is the same on each quadrature point). + * Stores whether a face (and both cells adjacent to the face) is + * Cartesian (face type 0), whether it represents an affine situation + * (face type 1), whether it is a flat face where the normal vector is + * the same throughout the face (face type 2), or is general (face type + * 3). */ - AlignedVector >, - VectorizedArray > > cartesian_data; + std::vector face_type; /** - * The first field stores the Jacobian for non-Cartesian cells where all - * the Jacobians on the cell are the same (i.e., constant, which comes - * from a linear transformation from unit to real cell). Also use - * similarities from one cell to another (on structured meshes, there - * are usually many cells with the same Jacobian). - * - * The second field stores the Jacobian determinant for non-Cartesian - * cells with constant Jacobian throughout the cell (without the - * quadrature weight, which depends on the quadrature point, whereas the - * determinant is the same on each quadrature point). + * The data cache for the cells. */ - AlignedVector >, - VectorizedArray > > affine_data; + std::vector > cell_data; /** - * Definition of a structure that stores data that depends on the - * quadrature formula (if we have more than one quadrature formula on a - * given problem, these fields will be different) + * The data cache for the faces. */ - struct MappingInfoDependent - { - /** - * This field stores the row starts for the inverse Jacobian - * transformations, quadrature weights and second derivatives. - */ - std::vector rowstart_jacobians; - - /** - * This field stores the inverse Jacobian transformation from unit to - * real cell, which is needed for most gradient transformations - * (corresponds to FEValues::inverse_jacobian) for general cells. - */ - AlignedVector > > jacobians; - - /** - * This field stores the Jacobian determinant times the quadrature - * weights (JxW in deal.II speak) for general cells. - */ - AlignedVector > JxW_values; - - /** - * Stores the diagonal part of the gradient of the inverse Jacobian - * transformation. The first index runs over the derivatives - * $\partial^2/\partial x_i^2$, the second over the space coordinate. - * Needed for computing the Laplacian of FE functions on the real - * cell. Uses a separate storage from the off-diagonal part - * $\partial^2/\partial x_i \partial x_j, i\neq j$ because that is - * only needed for computing a full Hessian. - */ - AlignedVector > > jacobians_grad_diag; - - /** - * Stores the off-diagonal part of the gradient of the inverse - * Jacobian transformation. Because of symmetry, only the upper - * diagonal part is needed. The first index runs through the - * derivatives row-wise, i.e., $\partial^2/\partial x_1 \partial x_2$ - * first, then $\partial^2/\partial x_1 \partial x_3$, and so on. The - * second index is the spatial coordinate. Not filled currently. - */ - AlignedVector1?dim*(dim-1)/2:1), - Tensor<1,dim,VectorizedArray > > > jacobians_grad_upper; + std::vector > face_data; - /** - * Stores the row start for quadrature points in real coordinates for - * both types of cells. Note that Cartesian cells will have shorter - * fields (length is @p n_q_points_1d) than non-Cartesian cells - * (length is @p n_q_points). - */ - std::vector rowstart_q_points; - - /** - * Stores the quadrature points in real coordinates for Cartesian - * cells (does not need to store the full data on all points) - */ - AlignedVector > > quadrature_points; - - /** - * The dim-dimensional quadrature formula underlying the problem - * (constructed from a 1D tensor product quadrature formula). - */ - dealii::hp::QCollection quadrature; + /** + * The data cache for the face-associated-with-cell topology, following + * the @p cell_type variable for the cell types. + */ + std::vector > face_data_by_cells; - /** - * The (dim-1)-dimensional quadrature formula corresponding to face - * evaluation (constructed from a 1D tensor product quadrature - * formula). - */ - dealii::hp::QCollection face_quadrature; + /** + * Computes the information in the given cells, called within + * initialize. + */ + void initialize_cells (const dealii::Triangulation &tria, + const std::vector > &cells, + const std::vector &active_fe_index, + const Mapping &mapping, + const std::vector > &quad, + const UpdateFlags update_flags_cells); - /** - * The number of quadrature points for the current quadrature formula. - */ - std::vector n_q_points; + /** + * Helper function to determine which update flags must be set in the + * internal functions to initialize all data as requested by the user. + */ + static UpdateFlags + compute_update_flags (const UpdateFlags update_flags, + const std::vector > &quad = + std::vector >()); + }; - /** - * The number of quadrature points for the current quadrature formula - * when applied to a face. Only set if the quadrature formula is - * derived from a tensor product, since it is not defined from the - * full quadrature formula otherwise. - */ - std::vector n_q_points_face; - /** - * The quadrature weights (vectorized data format) on the unit cell. - */ - std::vector > > quadrature_weights; - /** - * This variable stores the number of quadrature points for all - * quadrature indices in the underlying element for easier access to - * data in the hp case. - */ - std::vector quad_index_conversion; + /** + * A helper class to extract either cell or face data from mapping info + * for use in FEEvaluationBase. + * + * @author Katharina Kormann, Martin Kronbichler, 2018 + */ + template struct MappingInfoCellsOrFaces; - /** - * Return the quadrature index for a given number of quadrature - * points. If not in hp mode or if the index is not found, this - * function always returns index 0. Hence, this function does not - * check whether the given degree is actually present. - */ - unsigned int - quad_index_from_n_q_points (const unsigned int n_q_points) const; + template + struct MappingInfoCellsOrFaces + { + static const MappingInfoStorage * + get(const MappingInfo &mapping_info, + const unsigned int quad_no) + { + AssertIndexRange(quad_no, mapping_info.cell_data.size()); + return &mapping_info.cell_data[quad_no]; + } + }; + template + struct MappingInfoCellsOrFaces + { + static const MappingInfoStorage * + get(const MappingInfo &mapping_info, + const unsigned int quad_no) + { + AssertIndexRange(quad_no, mapping_info.face_data.size()); + return &mapping_info.face_data[quad_no]; + } + }; - /** - * Prints a detailed summary of memory consumption in the different - * structures of this class to the given output stream. - */ - template - void print_memory_consumption(StreamType &out, - const SizeInfo &size_info) const; - /** - * Return the memory consumption in bytes. - */ - std::size_t memory_consumption () const; - }; - /** - * Contains all the stuff that depends on the quadrature formula - */ - std::vector mapping_data_gen; + /** + * A class that is used to compare floating point arrays (e.g. std::vectors, + * Tensor<1,dim>, etc.). The idea of this class is to consider two arrays as + * equal if they are the same within a given tolerance. We use this + * comparator class within an 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 + struct FPArrayComparator + { + FPArrayComparator (const Number scaling); - /** - * Stores whether JxW values have been initialized - */ - bool JxW_values_initialized; + bool operator() (const std::vector &v1, + const std::vector &v2) const; - /** - * Stores whether we computed second derivatives. - */ - bool second_derivatives_initialized; + bool operator ()(const Tensor<1,VectorizedArray::n_array_elements,Number> &t1, + const Tensor<1,VectorizedArray::n_array_elements,Number> &t2) const; - /** - * Stores whether we computed quadrature points. - */ - bool quadrature_points_initialized; + template + bool operator ()(const Tensor<1,dim,Tensor<1,VectorizedArray::n_array_elements,Number> > &t1, + const Tensor<1,dim,Tensor<1,VectorizedArray::n_array_elements,Number> > &t2) const; - /** - * Internal temporary data used for the initialization. - */ - struct CellData - { - CellData (const double jac_size); - void resize (const unsigned int size); - - AlignedVector > > quadrature_points; - AlignedVector > > general_jac; - AlignedVector > > general_jac_grad; - Tensor<2,dim,VectorizedArray > const_jac; - const double jac_size; - }; + template + bool operator ()(const Tensor<2,dim,Tensor<1,VectorizedArray::n_array_elements,Number> > &t1, + const Tensor<2,dim,Tensor<1,VectorizedArray::n_array_elements,Number> > &t2) const; - /** - * Helper function called internally during the initialize function. - */ - void evaluate_on_cell (const dealii::Triangulation &tria, - const std::pair *cells, - const unsigned int cell, - const unsigned int my_q, - CellType (&cell_t_prev)[n_vector_elements], - CellType (&cell_t)[n_vector_elements], - dealii::FEValues &fe_values, - CellData &cell_data) const; + Number tolerance; }; /* ------------------- inline functions ----------------------------- */ - template + template inline unsigned int - MappingInfo::MappingInfoDependent:: - quad_index_from_n_q_points (const unsigned int n_q_points) const + MappingInfoStorage + ::quad_index_from_n_q_points (const unsigned int n_q_points) const { - for (unsigned int i=0; i::get_cell_type (const unsigned int cell_no) const { AssertIndexRange (cell_no, cell_type.size()); - CellType enum_cell_type = (CellType)(cell_type[cell_no] % n_cell_types); - Assert(enum_cell_type != undefined, ExcInternalError()); - return enum_cell_type; - } - - - - template - inline - unsigned int - MappingInfo::get_cell_data_index (const unsigned int cell_no) const - { - AssertIndexRange (cell_no, cell_type.size()); - return cell_type[cell_no] >> n_cell_type_bits; + return cell_type[cell_no]; } } // end of namespace MatrixFreeFunctions diff --git a/include/deal.II/matrix_free/mapping_info.templates.h b/include/deal.II/matrix_free/mapping_info.templates.h index 8c79097a3f..f08e728baf 100644 --- a/include/deal.II/matrix_free/mapping_info.templates.h +++ b/include/deal.II/matrix_free/mapping_info.templates.h @@ -18,6 +18,8 @@ #include #include +#include +#include #include #include #include @@ -32,14 +34,165 @@ namespace internal { namespace MatrixFreeFunctions { - // ----------------- actual MappingInfo functions ------------------------- + + /* ------------------------ MappingInfoStorage implementation ---------- */ + + template + MappingInfoStorage::QuadratureDescriptor + ::QuadratureDescriptor() + : + n_q_points (numbers::invalid_unsigned_int) + { + } + + + + template + void + MappingInfoStorage::QuadratureDescriptor + ::initialize(const Quadrature<1> &quadrature_1d, + const UpdateFlags update_flags_inner_faces) + { + Assert(structdim+1 <= spacedim || + update_flags_inner_faces == update_default, + ExcMessage("Volume cells do not allow for setting inner faces")); + quadrature = Quadrature(quadrature_1d); + n_q_points = quadrature.size(); + quadrature_weights.resize(n_q_points); + for (unsigned int i=0; i + std::size_t + MappingInfoStorage::QuadratureDescriptor + ::memory_consumption() const + { + std::size_t memory = + sizeof (this) + + quadrature.memory_consumption() + + quadrature_weights.memory_consumption() + + face_orientations.memory_consumption(); + for (unsigned int d=0; d + std::size_t + MappingInfoStorage::memory_consumption() const + { + return + MemoryConsumption::memory_consumption (descriptor) + + MemoryConsumption::memory_consumption (data_index_offsets) + + MemoryConsumption::memory_consumption (JxW_values) + + MemoryConsumption::memory_consumption (normal_vectors) + + MemoryConsumption::memory_consumption (jacobians[0]) + + MemoryConsumption::memory_consumption (jacobians[1]) + + MemoryConsumption::memory_consumption (jacobian_gradients[0]) + + MemoryConsumption::memory_consumption (jacobian_gradients[1]) + + MemoryConsumption::memory_consumption (normals_times_jacobians[0]) + + MemoryConsumption::memory_consumption (normals_times_jacobians[1]) + + MemoryConsumption::memory_consumption (quadrature_point_offsets) + + MemoryConsumption::memory_consumption (quadrature_points); + } + + + + template + template + void + MappingInfoStorage::print_memory_consumption + (StreamType &out, + const SizeInfo &task_info) const + { + // print_memory_statistics involves global communication, so we can + // disable the check here only if no processor has any such data + const std::size_t size = Utilities::MPI::sum(jacobians[0].size(), + task_info.communicator); + if (size > 0) + { + out << " Memory JxW data: "; + task_info.print_memory_statistics + (out, MemoryConsumption::memory_consumption (data_index_offsets) + + MemoryConsumption::memory_consumption (JxW_values)); + out << " Memory Jacobian data: "; + task_info.print_memory_statistics + (out, MemoryConsumption::memory_consumption (jacobians[0]) + + MemoryConsumption::memory_consumption (jacobians[1])); + out << " Memory second derivative data: "; + task_info.print_memory_statistics + (out, MemoryConsumption::memory_consumption (jacobian_gradients[0]) + + MemoryConsumption::memory_consumption (jacobian_gradients[1])); + } + const std::size_t normal_size = Utilities::MPI::sum(normal_vectors.size(), + task_info.communicator); + if (normal_size > 0) + { + out << " Memory normal vectors data: "; + task_info.print_memory_statistics + (out, MemoryConsumption::memory_consumption (normal_vectors) + + MemoryConsumption::memory_consumption (normals_times_jacobians[0]) + + MemoryConsumption::memory_consumption (normals_times_jacobians[1])); + } + + const std::size_t quad_size = + Utilities::MPI::sum(quadrature_points.size(), + task_info.communicator); + if (quad_size > 0) + { + out << " Memory quadrature points: "; + task_info.print_memory_statistics + (out, MemoryConsumption::memory_consumption (quadrature_point_offsets) + + MemoryConsumption::memory_consumption (quadrature_points)); + } + } + + + + /* ------------------------ MappingInfo implementation ----------------- */ template MappingInfo::MappingInfo() - : - JxW_values_initialized (false), - second_derivatives_initialized (false), - quadrature_points_initialized (false) {} @@ -48,13 +201,11 @@ namespace internal void MappingInfo::clear () { - JxW_values_initialized = false; - quadrature_points_initialized = false; - second_derivatives_initialized = false; - mapping_data_gen.clear(); + cell_data.clear(); + face_data.clear(); + face_data_by_cells.clear(); cell_type.clear(); - cartesian_data.clear(); - affine_data.clear(); + face_type.clear(); } @@ -65,17 +216,10 @@ namespace internal compute_update_flags (const UpdateFlags update_flags, const std::vector > &quad) { - // this class is build around the evaluation this class is build around - // the evaluation of inverse gradients, so compute them in any case - UpdateFlags new_flags = update_inverse_jacobians; - - // if the user requested gradients, need inverse Jacobians - if (update_flags & update_gradients || update_flags & update_inverse_jacobians) - new_flags |= update_inverse_jacobians; - - // for JxW, would only need JxW values. - if (update_flags & update_JxW_values) - new_flags |= update_JxW_values; + // this class is build around the evaluation of jacobians, so compute + // them in any case. The Jacobians will be inverted manually. Since we + // always do support integration, we also include the JxW values + UpdateFlags new_flags = update_jacobians | update_JxW_values; // for Hessian information, need inverse Jacobians and the derivative of // Jacobians (these two together will give use the gradients of the @@ -112,7 +256,30 @@ namespace internal - namespace internal + template + void + MappingInfo::initialize + (const dealii::Triangulation &tria, + const std::vector > &cells, + const std::vector &active_fe_index, + const Mapping &mapping, + const std::vector > &quad, + const UpdateFlags update_flags_cells) + { + clear(); + + // 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, quad, update_flags_cells); + } + + + + /* ------------------------- initialization of cells ------------------- */ + + // Anonymous namespace with implementation of extraction of values on cell + // range + namespace { template double get_jacobian_size (const dealii::Triangulation &tria) @@ -121,331 +288,438 @@ namespace internal return 1; else return tria.begin()->diameter(); } - } - template - void - MappingInfo::initialize - (const dealii::Triangulation &tria, - const std::vector > &cells, - const std::vector &active_fe_index, - const Mapping &mapping, - const std::vector > &quad, - const UpdateFlags update_flags_input) - { - clear(); - const unsigned int n_quads = quad.size(); - const unsigned int n_cells = cells.size(); - const unsigned int vectorization_length = - VectorizedArray::n_array_elements; - Assert (n_cells%vectorization_length == 0, ExcInternalError()); - const unsigned int n_macro_cells = n_cells/vectorization_length; - mapping_data_gen.resize (n_quads); - cell_type.resize (n_macro_cells); + template + struct CompressedCellData + { + CompressedCellData(const double expected_size) + : + data(FPArrayComparator(expected_size)) + {} + + std::map::n_array_elements,Number> > + , unsigned int, FPArrayComparator > data; + }; + + /** + * Internal temporary data used for the initialization. + */ + template + struct LocalData + { + LocalData (const double jac_size); + void resize (const unsigned int size); - // dummy FE that is used to set up an FEValues object. Do not need the - // actual finite element because we will only evaluate quantities for - // the mapping that are independent of the FE - FE_Nothing dummy_fe; - UpdateFlags update_flags = compute_update_flags (update_flags_input, quad); + AlignedVector > > quadrature_points; + AlignedVector > > general_jac; + AlignedVector > JxW_values; + AlignedVector > > general_jac_grad; + AlignedVector > > normal_vectors; + Tensor<2,dim,VectorizedArray > const_jac; + const double jac_size; + }; - if (update_flags & update_JxW_values) - JxW_values_initialized = true; - if (update_flags & update_jacobian_grads) - second_derivatives_initialized = true; - if (update_flags & update_quadrature_points) - quadrature_points_initialized = true; - - // 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 of the - // first cell (on level zero, which is the best accuracy we can hope - // for, since diameters on finer levels are computed by differences of - // nearby cells). If the mesh extends over a certain domain, the - // precision of double values is essentially limited by this precision. - const double jacobian_size = internal::get_jacobian_size(tria); - - // objects that hold the data for up to vectorization_length cells while - // we fill them up. Only after all vectorization_length cells have been - // processed, we can insert the data into the data structures of this - // class - CellData data (jacobian_size); - for (unsigned int my_q=0; my_q n_q_points_1d (n_hp_quads), - step_size_cartesian (n_hp_quads); - if (n_hp_quads > 1) - current_data.quad_index_conversion.resize(n_hp_quads); - for (unsigned int q=0; q(n_q_points_1d[q]); - current_data.n_q_points.push_back (n_q_points); - - current_data.n_q_points_face.push_back - (dim>1 ? Utilities::fixed_power(n_q_points_1d[q]) : 1); - current_data.quadrature.push_back - (Quadrature(quad[my_q][q])); - current_data.face_quadrature.push_back - (Quadrature(quad[my_q][q])); - - // set quadrature weights in vectorized form - current_data.quadrature_weights[q].resize(n_q_points); - for (unsigned int i=0; i 1) - current_data.quad_index_conversion[q] = n_q_points; - - // 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[q] = 0; - unsigned int factor = 1; - for (unsigned int d=0; d + LocalData::LocalData (const double jac_size_in) + : + jac_size (jac_size_in) + {} + + + + template + void + LocalData::resize (const unsigned int size) + { + if (JxW_values.size() != size) + { + quadrature_points.resize_fast(size); + general_jac.resize_fast(size*2); + JxW_values.resize_fast(size); + general_jac_grad.resize_fast(size*2); + normal_vectors.resize_fast(size); + } + } + + /** + * Helper function called internally during the initialize function. + */ + template + void evaluate_on_cell (const dealii::Triangulation &tria, + const std::pair *cells, + const unsigned int my_q, + CellType &cell_t_prev, + CellType (&cell_t)[VectorizedArray::n_array_elements], + dealii::FEValues &fe_val, + LocalData &cell_data) + { + const unsigned int n_q_points = fe_val.n_quadrature_points; + const UpdateFlags update_flags = fe_val.get_update_flags(); + + cell_data.const_jac = Tensor<2,dim,VectorizedArray >(); + + // this should be the same value as used in HashValue::scaling (but we + // not have that field here) + const double zero_tolerance_double = cell_data.jac_size * + std::numeric_limits::epsilon() * 1024.; + for (unsigned int j=0; j::n_array_elements; ++j) + { + typename dealii::Triangulation::cell_iterator + cell_it (&tria, cells[j].first, cells[j].second); + fe_val.reinit(cell_it); + cell_t[j] = general; + + // extract quadrature points and store them temporarily. if we have + // Cartesian cells, we can compress the indices + if (update_flags & update_quadrature_points) + for (unsigned int q=0; q &point = fe_val.quadrature_point(q); + for (unsigned int d=0; d > jac_grad, grad_jac_inv; - Tensor<1,dim,VectorizedArray > tmp; - - // encodes the cell types of the current cell. Since several cells - // must be considered together, this variable holds the individual - // info of the last chunk of cells - CellType cell_t [vectorization_length], - cell_t_prev [vectorization_length]; - for (unsigned int j=0; j::active_cell_iterator, we need to manually - // select the correct finite element, so just hold a vector of - // FEValues - std::vector > > - fe_values (current_data.quadrature.size()); - UpdateFlags update_flags_feval = - ((update_flags & update_inverse_jacobians) ? update_jacobians : update_default) | - ((update_flags & update_jacobian_grads) ? update_jacobian_grads : update_default) | - ((update_flags & update_quadrature_points) ? update_quadrature_points : update_default); - - // resize the fields that have fixed size or for which we know - // something from an earlier loop - current_data.rowstart_q_points.resize (n_macro_cells+1); - if (my_q > 0) - { - const unsigned int n_cells_var = - mapping_data_gen[0].rowstart_jacobians.size()-1; - current_data.rowstart_jacobians.reserve (n_cells_var+1); - const unsigned int reserve_size = n_cells_var * - current_data.n_q_points[0]; - if (mapping_data_gen[0].jacobians.size() > 0) - current_data.jacobians.reserve (reserve_size); - if (mapping_data_gen[0].JxW_values.size() > 0) - current_data.jacobians.reserve (reserve_size); - if (mapping_data_gen[0].jacobians_grad_diag.size() > 0) - current_data.jacobians_grad_diag.reserve (reserve_size); - if (mapping_data_gen[0].jacobians_grad_upper.size() > 0) - current_data.jacobians_grad_upper.reserve (reserve_size); - } + // if this is not the first quadrature formula and we already have + // determined that this cell is either Cartesian or with constant + // Jacobian, we have nothing more to do. + if (my_q > 0 && cell_t_prev <= affine) + continue; + + // first round: if the transformation is detected to be the same as + // on the old cell, we only need to copy over the data. + if (fe_val.get_cell_similarity() == CellSimilarity::translation + && + my_q == 0) + { + if (j==0) + cell_t[j] = cell_t_prev; + else + cell_t[j] = cell_t[j-1]; + } + + const DerivativeForm<1,dim,dim> &jac_0 = fe_val.jacobian(0); + + if (my_q == 0) + { + // check whether the Jacobian is constant on this cell the first + // time we come around here + if (cell_t[j] == general) + { + bool jacobian_constant = true; + for (unsigned int q=1; q &jac = fe_val.jacobian(q); + for (unsigned int d=0; d + zero_tolerance_double) + jacobian_constant = false; + if (jacobian_constant == false) + break; + } + + // check whether the Jacobian is diagonal to machine + // accuracy + bool cell_cartesian = jacobian_constant; + for (unsigned int d=0; d + zero_tolerance_double) + { + cell_cartesian=false; + break; + } + + // in case we have only one quadrature point, we can have + // non-constant Jacobians, but we cannot detect it by + // comparison from one quadrature point to the next: in that + // case, need to look at second derivatives and see whether + // there are some non-zero entries (this is necessary since + // we determine the constness of the Jacobian for the first + // quadrature formula and might not look at them any more + // for the second, third quadrature formula). in any case, + // the flag update_jacobian_grads will be set in that case + if (cell_cartesian == false && n_q_points == 1 && + update_flags & update_jacobian_grads) + { + const DerivativeForm<1,dim,dim> &jac = fe_val.jacobian(0); + const DerivativeForm<2,dim,dim> &jacobian_grad = + fe_val.jacobian_grad(0); + for (unsigned int d=0; d + zero_tolerance_double) + jacobian_constant = false; + } + } + // set cell type + if (cell_cartesian == true) + cell_t[j] = cartesian; + else if (jacobian_constant == true) + cell_t[j] = affine; + else + cell_t[j] = general; + } + + // Cartesian cell + if (cell_t[j] == cartesian) + { + // set Jacobian into diagonal (off-diagonal part is already + // zeroed out) + for (unsigned int d=0; d > as - // key into the std::map, but std::map allocation does not align the - // allocated memory correctly, so put it into a tensor of the - // correct length instead - FPArrayComparator comparator(jacobian_size); - typedef Tensor<1,VectorizedArray::n_array_elements,Number> VEC_ARRAY; - std::map, unsigned int, - FPArrayComparator > cartesians(comparator); - std::map, unsigned int, - FPArrayComparator > affines(comparator); - - // loop over all cells - for (unsigned int cell=0; cell &jac = fe_val.jacobian(q); + for (unsigned int d=0; d &jacobian_grad = fe_val.jacobian_grad(q); + for (unsigned int d=0; d::n_array_elements-1]; + } + + + + template + void + initialize_cell_range + (const std::pair cell_range, + const dealii::Triangulation &tria, + const std::vector > &cells, + const std::vector &active_fe_index, + const Mapping &mapping, + const std::vector > &quad, + const UpdateFlags update_flags, + MappingInfo &mapping_info, + std::pair >, + CompressedCellData > &data) + { + FE_Nothing dummy_fe; + + Tensor<3,dim,VectorizedArray > jac_grad, grad_jac_inv; + Tensor<1,dim,VectorizedArray > 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 + // of the first cell (on level zero, which is the best accuracy we can + // hope for, since diameters on finer levels are computed by + // differences of nearby cells) as the order of magnitude by which we + // make comparisons "relative." + const double jacobian_size = get_jacobian_size(tria); + + // objects that hold the data for up to vectorization_width cells while + // we fill them up. Only after all vectorization_width cells have been + // processed, we can insert the data into the data structures of this + // class + LocalData cell_data (jacobian_size); + + // encodes the cell types of the current cell. Since several cells + // must be considered together, this variable holds the individual + // info of the last chunk of cells + CellType cell_t [VectorizedArray::n_array_elements]; + CellType cell_t_prev = general; + + // fe_values object that is used to compute the mapping data. for + // the hp case there might be more than one finite element. since we + // manually select the active FE index and not via a + // hp::DoFHandler::active_cell_iterator, we need to manually + // select the correct finite element, so just hold a vector of + // FEValues + std::vector > > > + fe_values (mapping_info.cell_data.size()); + for (unsigned int i=0; i > n_q_points_1d (quad.size()), + step_size_cartesian (quad.size()); + for (unsigned int my_q=0; my_q 0 ? active_fe_index[cell] : 0; - const unsigned int n_q_points = current_data.n_q_points[fe_index]; - if (fe_values[fe_index].get() == nullptr) - fe_values[fe_index].reset - (new dealii::FEValues (mapping, dummy_fe, - current_data.quadrature[fe_index], + const unsigned int n_q_points = + mapping_info.cell_data[my_q].descriptor[fe_index].n_q_points; + if (fe_values[my_q][fe_index].get() == nullptr) + fe_values[my_q][fe_index].reset + (new dealii::FEValues (mapping, dummy_fe, mapping_info.cell_data[my_q]. + descriptor[fe_index].quadrature, update_flags_feval)); - dealii::FEValues &fe_val = *fe_values[fe_index]; - data.resize (n_q_points); + dealii::FEValues &fe_val = *fe_values[my_q][fe_index]; + cell_data.resize (n_q_points); // if the fe index has changed from the previous cell, set the // old cell type to invalid (otherwise, we might detect // similarity due to some cells further ahead) - if (cell > 0 && active_fe_index.size() > 0 && - active_fe_index[cell] != active_fe_index[cell-1]) - cell_t_prev[vectorization_length-1] = undefined; - evaluate_on_cell (tria, &cells[cell*vectorization_length], - cell, my_q, cell_t_prev, cell_t, fe_val, data); + if (my_q > 0) + cell_t_prev = CellType(mapping_info.cell_type[cell]); + else if (cell > cell_range.first && active_fe_index.size() > 0 && + active_fe_index[cell] != active_fe_index[cell-1]) + cell_t_prev = general; + + evaluate_on_cell (tria, &cells[cell*VectorizedArray::n_array_elements], + my_q, cell_t_prev, cell_t, fe_val, + cell_data); // now reorder the data into vectorized types. if we are here // for the first time, we need to find out whether the Jacobian // allows for some simplification (Cartesian, affine) taking - // vectorization_length cell together and we have to insert that - // data into the respective fields. Also, we have to compress - // different cell indicators into one structure. + // vectorization_width cell together if (my_q == 0) { - // find the most general cell type (most general type is 2 + // find the most general cell type (most general type is 3 // (general cell)) CellType most_general_type = cartesian; - for (unsigned int j=0; j::n_array_elements; ++j) if (cell_t[j] > most_general_type) most_general_type = cell_t[j]; - AssertIndexRange (most_general_type, 3); - unsigned int insert_position = numbers::invalid_unsigned_int; - - // Cartesian cell with diagonal Jacobian: only insert the - // diagonal of the inverse and the Jacobian determinant. We - // do this by using an std::map that collects pointers to - // all Cartesian Jacobians. We need a pointer in the - // std::map because it cannot store data based on - // VectorizedArray (alignment issue). We circumvent the - // problem by temporarily filling the next value into the - // cartesian_data field and, in case we did an insertion, - // the data is already in the correct place. - if (most_general_type == cartesian) - { - std::pair,unsigned int> new_entry; - new_entry.second = cartesians.size(); - for (unsigned int d=0; d::n_array_elements; ++v) - new_entry.first[d][v] = data.const_jac[d][d][v]; - - std::pair, - unsigned int, FPArrayComparator >::iterator, - bool> it = cartesians.insert(new_entry); - if (it.second == false) - insert_position = it.first->second; - else - insert_position = new_entry.second; - } + AssertIndexRange ((unsigned int)most_general_type, 4U); + mapping_info.cell_type[cell] = most_general_type; + } + + AssertThrow(data.first[my_q].JxW_values.size() < + static_cast(std::numeric_limits::max()), + ExcMessage("Index overflow. Cannot fit data in 32 bit integers")); - // Constant Jacobian case. same strategy as before, but with - // other data fields - else if (most_general_type == affine) + unsigned int insert_position = data.first[my_q].JxW_values.size(); + // Cartesian/affine cell with constant Jacobians throughout the + // cell. We need to store the data in another data field because + // std::map cannot store data based on VectorizedArray directly + // (alignment issue). + if (mapping_info.cell_type[cell] <= affine) + { + if (my_q == 0) { - std::pair,unsigned int> new_entry; - new_entry.second = affines.size(); + std::pair:: + n_array_elements,Number> >,unsigned int> new_entry; + // This number overlaps with the general data but we + // take care of that when we merge data from different + // threads + new_entry.second = data.second.data.size(); for (unsigned int d=0; d::n_array_elements; ++v) - new_entry.first[d][e][v] = data.const_jac[d][e][v]; - - std::pair, - unsigned int, FPArrayComparator >::iterator, - bool> it = affines.insert(new_entry); - if (it.second == false) - insert_position = it.first->second; - else - insert_position = new_entry.second; - } + new_entry.first[d][e][v] = cell_data.const_jac[d][e][v]; - // general cell case: first resize the data field to fit the - // new data. if we are here the first time, assume that - // there are many general cells to come, so reserve some - // memory in order to not have too many reallocations and - // memcpy's. The scheme used here involves at most one - // reallocation. - else - { - Assert (most_general_type == general, ExcInternalError()); - insert_position = current_data.rowstart_jacobians.size(); - if (current_data.rowstart_jacobians.size() == 0) - { - unsigned int reserve_size = (n_macro_cells-cell+1)/2; - current_data.rowstart_jacobians.reserve - (reserve_size); - reserve_size *= n_q_points; - current_data.jacobians.reserve (reserve_size); - if (update_flags & update_JxW_values) - current_data.JxW_values.reserve (reserve_size); - if (update_flags & update_jacobian_grads) - { - current_data.jacobians_grad_diag.reserve (reserve_size); - current_data.jacobians_grad_upper.reserve (reserve_size); - } - } + insert_position = data.second.data.insert(new_entry).first->second; } - - cell_type[cell] = ((insert_position << n_cell_type_bits) + - (unsigned int)most_general_type); - - } // end if (my_q == 0) + else + insert_position = data.first[0].data_index_offsets[cell-cell_range.first]; + } // general cell case: now go through all quadrature points and // collect the data. done for all different quadrature formulas, // so do it outside the above loop. - if (get_cell_type(cell) == general) + data.first[my_q].data_index_offsets.push_back(insert_position); + if (mapping_info.get_cell_type(cell) == general) { - const unsigned int previous_size = - current_data.jacobians.size(); - current_data.rowstart_jacobians.push_back (previous_size); - if (update_flags & update_JxW_values) - { - AssertDimension (previous_size, - current_data.JxW_values.size()); - } - if (update_flags & update_jacobian_grads) - { - AssertDimension (previous_size, - current_data.jacobians_grad_diag.size()); - AssertDimension (previous_size, - current_data.jacobians_grad_upper.size()); - } for (unsigned int q=0; q > &jac = data.general_jac[q]; - Tensor<3,dim,VectorizedArray > &jacobian_grad = data.general_jac_grad[q]; - for (unsigned int j=0; j > &jac = cell_data.general_jac[q]; + Tensor<3,dim,VectorizedArray > &jacobian_grad = cell_data.general_jac_grad[q]; + for (unsigned int j=0; j::n_array_elements; ++j) + if (cell_t[j] < general) { for (unsigned int d=0; d det = determinant (jac); - current_data.jacobians.push_back (transpose(invert(jac))); - const Tensor<2,dim,VectorizedArray > &inv_jac = current_data.jacobians.back(); - - if (update_flags & update_JxW_values) - current_data.JxW_values.push_back - (det * current_data.quadrature_weights[fe_index][q]); + data.first[my_q].JxW_values.push_back(determinant(jac)* + fe_val.get_quadrature().weight(q)); + Tensor<2,dim,VectorizedArray > inv_jac = transpose(invert(jac)); + data.first[my_q].jacobians[0].push_back(inv_jac); if (update_flags & update_jacobian_grads) { @@ -498,22 +772,18 @@ namespace internal } } - { - VectorizedArray grad_diag[dim][dim]; - for (unsigned int d=0; d >(grad_diag)); - } + // the diagonal part of Jacobian gradient comes first + Tensor<1,dim*(dim+1)/2,Tensor<1,dim,VectorizedArray > > final_grad; + for (unsigned int d=0; d1?dim*(dim-1)/2:1),Tensor<1,dim,VectorizedArray > > grad_upper; + // then the upper-diagonal part for (unsigned int d=0, count=0; d > quad_point; - - if (get_cell_type(cell) == cartesian) + if (mapping_info.get_cell_type(cell) == cartesian) { - current_data.quadrature_points.resize (old_size+ - n_q_points_1d[fe_index]); - for (unsigned int q=0; q > quad_point; + for (unsigned int d=0; d, - unsigned int, FPArrayComparator >::iterator - it = cartesians.begin(); it != cartesians.end(); ++it) - { - VectorizedArray det = make_vectorized_array(1.); - for (unsigned int d=0; d jac_d; - for (unsigned int v=0; - v::n_array_elements; ++v) - jac_d[v] = it->first[d][v]; - cartesian_data[it->second].first[d] = 1./jac_d; - det *= jac_d; - } - cartesian_data[it->second].second = det; - } - affine_data.resize(affines.size()); - for (typename std::map, - unsigned int, FPArrayComparator >::iterator - it = affines.begin(); it != affines.end(); ++it) - { - Tensor<2,dim,VectorizedArray > jac; - for (unsigned int d=0; d::n_array_elements; ++v) - jac[d][e][v] = it->first[d][e][v]; + } // end for ( cell < end_cells ) + } - affine_data[it->second].second = determinant(jac); - affine_data[it->second].first = transpose(invert(jac)); - } - } - } - } + + + template + void + merge_compressed_data(const CONTAINER &source, + CONTAINER &destination, + std::vector &indices) + { + indices.resize(source.size()); + typename CONTAINER::iterator lookup = destination.begin(); + for (typename CONTAINER::const_iterator it = source.begin(); + it != source.end(); ++it) + { + typename CONTAINER::value_type entry = *it; + entry.second = destination.size(); + lookup = destination.insert(lookup, entry); + AssertIndexRange(it->second, indices.size()); + indices[it->second] = lookup->second; + // best guess for insert position of next item + ++lookup; + } + } + + + + template + void + copy_data (const unsigned int first_cell, + const std::array &data_shift, + const std::vector &indices_compressed, + const std::vector &cell_type, + MappingInfoStorage &data_cells_local, + MappingInfoStorage &data_cells) + { + // Copy the index offsets and shift by the appropriate value + for (unsigned int lcell=0; + lcell(affine) ? + 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] = + data_cells_local.quadrature_point_offsets[lcell] + + data_shift[1]; + } + + // Copy quadrature points + if (data_cells.quadrature_point_offsets.empty() == false) + { + Point > *out_point = + &data_cells.quadrature_points[data_shift[1]]; + for (const Point > *point = + data_cells_local.quadrature_points.begin(); point != + data_cells_local.quadrature_points.end(); ++point, ++out_point) + *out_point = *point; + data_cells_local.quadrature_points.clear(); + } + + // If we have collected Jacobian data, copy Jacobians, JxW values, + // Jacobian gradients + if (data_cells_local.JxW_values.empty()) + return; + + std::copy(data_cells_local.JxW_values.begin(), + data_cells_local.JxW_values.end(), + data_cells.JxW_values.begin()+data_shift[0]); + data_cells_local.JxW_values.clear(); + std::copy(data_cells_local.normal_vectors.begin(), + data_cells_local.normal_vectors.end(), + data_cells.normal_vectors.begin()+data_shift[0]); + data_cells_local.normal_vectors.clear(); + for (unsigned int i=0; i<2; ++i) + { + std::copy(data_cells_local.jacobians[i].begin(), + data_cells_local.jacobians[i].end(), + data_cells.jacobians[i].begin()+data_shift[0]); + data_cells_local.jacobians[i].clear(); + std::copy(data_cells_local.jacobian_gradients[i].begin(), + data_cells_local.jacobian_gradients[i].end(), + data_cells.jacobian_gradients[i].begin()+data_shift[0]); + data_cells_local.jacobian_gradients[i].clear(); + std::copy(data_cells_local.normals_times_jacobians[i].begin(), + data_cells_local.normals_times_jacobians[i].end(), + data_cells.normals_times_jacobians[i].begin()+data_shift[0]); + data_cells_local.normals_times_jacobians[i].clear(); + } + } + + } // end of anonymous namespace template void - MappingInfo::evaluate_on_cell (const dealii::Triangulation &tria, - const std::pair *cells, - const unsigned int cell, - const unsigned int my_q, - CellType (&cell_t_prev)[n_vector_elements], - CellType (&cell_t)[n_vector_elements], - dealii::FEValues &fe_val, - CellData &data) const + MappingInfo::initialize_cells + (const dealii::Triangulation &tria, + const std::vector > &cells, + const std::vector &active_fe_index, + const Mapping &mapping, + const std::vector > &quad, + const UpdateFlags update_flags_input) { - const unsigned int n_q_points = fe_val.n_quadrature_points; - const UpdateFlags update_flags = fe_val.get_update_flags(); - - // this should be the same value as used in HashValue::scaling (but we - // not have that field here) - const double zero_tolerance_double = data.jac_size * - std::numeric_limits::epsilon() * 1024.; - for (unsigned int j=0; j::cell_iterator - cell_it (&tria, cells[j].first, cells[j].second); - fe_val.reinit(cell_it); - cell_t[j] = undefined; + const unsigned int n_quads = quad.size(); + const unsigned int n_cells = cells.size(); + const unsigned int vectorization_width = + VectorizedArray::n_array_elements; + Assert (n_cells%vectorization_width == 0, ExcInternalError()); + const unsigned int n_macro_cells = n_cells/vectorization_width; + cell_data.resize (n_quads); + cell_type.resize (n_macro_cells); - // extract quadrature points and store them temporarily. if we have - // Cartesian cells, we can compress the indices - if (update_flags & update_quadrature_points) - for (unsigned int q=0; q &point = fe_val.quadrature_point(q); - for (unsigned int d=0; d 0 && (get_cell_type(cell) == cartesian - || get_cell_type(cell) == affine) ) - continue; - - // first round: if the transformation is detected to be the same as - // on the old cell, we only need to copy over the data. - if (fe_val.get_cell_similarity() == CellSimilarity::translation - && - my_q == 0) - { - if (j==0) - { - Assert (cell>0, ExcInternalError()); - cell_t[j] = cell_t_prev[n_vector_elements-1]; - } - else - cell_t[j] = cell_t[j-1]; - } + for (unsigned int my_q=0; my_q &jac_0 = fe_val.jacobian(0); + if (n_macro_cells == 0) + return; - if (my_q == 0) - { - // check whether the Jacobian is constant on this cell the first - // time we come around here - if (cell_t[j] == undefined) - { - bool jacobian_constant = true; - for (unsigned int q=1; q &jac = fe_val.jacobian(q); - for (unsigned int d=0; d - zero_tolerance_double) - jacobian_constant = false; - if (jacobian_constant == false) - break; - } + // Create as many chunks of cells as we have threads and spawn the work + unsigned int work_per_chunk = + std::max(8U, (n_macro_cells + MultithreadInfo::n_threads() - 1) / + MultithreadInfo::n_threads()); - // check whether the Jacobian is diagonal to machine - // accuracy - bool cell_cartesian = jacobian_constant; - for (unsigned int d=0; d - zero_tolerance_double) - { - cell_cartesian=false; - break; - } + std::vector >, + CompressedCellData > > data_cells_local; + // Reserve enough space to avoid re-allocation (which would break the + // references to the data fields passed to the tasks!) + data_cells_local.reserve(MultithreadInfo::n_threads()); - // in case we have only one quadrature point, we can have - // non-constant Jacobians, but we cannot detect it by - // comparison from one quadrature point to the next: in that - // case, need to look at second derivatives and see whether - // there are some non-zero entries (this is necessary since - // we determine the constness of the Jacobian for the first - // quadrature formula and might not look at them any more - // for the second, third quadrature formula). in any case, - // the flag update_jacobian_grads will be set in that case - if (cell_cartesian == false && n_q_points == 1 && - update_flags & update_jacobian_grads) - { - const DerivativeForm<1,dim,dim> &jac = fe_val.jacobian(0); - const DerivativeForm<2,dim,dim> &jacobian_grad = - fe_val.jacobian_grad(0); - for (unsigned int d=0; d - zero_tolerance_double) - jacobian_constant = false; - } - } - // set cell type - if (cell_cartesian == true) - cell_t[j] = cartesian; - else if (jacobian_constant == true) - cell_t[j] = affine; - else - cell_t[j] = general; - } + { + Threads::TaskGroup<> tasks; + std::pair cell_range(0U, work_per_chunk); + while (cell_range.first < n_macro_cells) + { + data_cells_local.push_back + (std::make_pair (std::vector >(n_quads), + CompressedCellData(get_jacobian_size(tria)))); + tasks += Threads::new_task(&initialize_cell_range, + cell_range, tria, + cells, active_fe_index, mapping, quad, + update_flags, *this, + data_cells_local.back()); + cell_range.first = cell_range.second; + cell_range.second += work_per_chunk; + } + tasks.join_all(); + } - // Cartesian cell - if (cell_t[j] == cartesian) - { - // set Jacobian into diagonal and clear off-diagonal part - for (unsigned int d=0; d > indices_compressed(data_cells_local.size()); + for (unsigned int i=0; i > shift(data_cells_local.size()); + shift[0][0] = n_constant_jacobians; + shift[0][1] = 0; + for (unsigned int i=1; i tasks; + for (unsigned int i=0; i, + work_per_chunk * i, shift[i], + indices_compressed[i], cell_type, + data_cells_local[i].first[my_q], + cell_data[my_q]); + + // finally, insert the constant cell data at the beginning (the + // other tasks can already start copying the non-constant data) + if (my_q == 0) + { + for (auto &it : data_cells_local[0].second.data) { - // compress out very small values + Tensor<2,dim,VectorizedArray > jac; for (unsigned int d=0; d::n_array_elements; ++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); + // invert and transpose jac + cell_data[my_q].jacobians[0][index] = transpose(invert(jac)); + // second derivative of transformation is zero on affine cells } } - - // general cell case - - // go through all quadrature points and fill in the data into the - // temporary data structures with slots for the vectorized data - // types - for (unsigned int q=0; q &jac = fe_val.jacobian(q); - for (unsigned int d=0; d &jacobian_grad = fe_val.jacobian_grad(q); - for (unsigned int d=0; d - MappingInfo::CellData::CellData (const double jac_size_in) - : - jac_size (jac_size_in) - {} + std::size_t MappingInfo::memory_consumption() const + { + std::size_t + memory = MemoryConsumption::memory_consumption (cell_data); + memory += MemoryConsumption::memory_consumption (face_data); + memory += cell_type.capacity()*sizeof(CellType); + memory += face_type.capacity()*sizeof(CellType); + memory += sizeof (*this); + return memory; + } template - void - MappingInfo::CellData::resize (const unsigned int size) + template + void MappingInfo::print_memory_consumption(StreamType &out, + const SizeInfo &task_info) const { - if (general_jac.size() != size) + out << " Cell types: "; + task_info.print_memory_statistics + (out, cell_type.capacity()*sizeof(CellType)); + out << " Face types: "; + task_info.print_memory_statistics + (out, face_type.capacity()*sizeof(CellType)); + for (unsigned int j=0; j - std::size_t MappingInfo::MappingInfoDependent::memory_consumption() const + /* ------------------------------------------------------------------ */ + + template + FPArrayComparator::FPArrayComparator (const Number scaling) + : + tolerance (scaling * std::numeric_limits::epsilon() * 1024.) + {} + + + + template + bool + FPArrayComparator::operator() (const std::vector &v1, + const std::vector &v2) const { - std::size_t - memory = MemoryConsumption::memory_consumption (jacobians); - memory += MemoryConsumption::memory_consumption (JxW_values); - memory += MemoryConsumption::memory_consumption (jacobians_grad_diag); - memory += MemoryConsumption::memory_consumption (jacobians_grad_upper); - memory += MemoryConsumption::memory_consumption (rowstart_q_points); - memory += MemoryConsumption::memory_consumption (quadrature_points); - memory += MemoryConsumption::memory_consumption (quadrature); - memory += MemoryConsumption::memory_consumption (face_quadrature); - memory += MemoryConsumption::memory_consumption (quadrature_weights); - memory += MemoryConsumption::memory_consumption (n_q_points); - memory += MemoryConsumption::memory_consumption (n_q_points_face); - memory += MemoryConsumption::memory_consumption (quad_index_conversion); - return memory; + const unsigned int s1 = v1.size(), s2 = v2.size(); + if (s1 < s2) + return true; + else if (s1 > s2) + return false; + else + for (unsigned int i=0; i v2[i] + tolerance) + return false; + return false; } - template - std::size_t MappingInfo::memory_consumption() const + template + bool + FPArrayComparator:: + operator ()(const Tensor<1,VectorizedArray::n_array_elements,Number> &t1, + const Tensor<1,VectorizedArray::n_array_elements,Number> &t2) const { - std::size_t - memory= MemoryConsumption::memory_consumption (mapping_data_gen); - memory += MemoryConsumption::memory_consumption (affine_data); - memory += MemoryConsumption::memory_consumption (cartesian_data); - memory += MemoryConsumption::memory_consumption (cell_type); - memory += sizeof (*this); - return memory; + for (unsigned int k=0; k::n_array_elements; ++k) + if (t1[k] < t2[k] - tolerance) + return true; + else if (t1[k] > t2[k] + tolerance) + return false; + return false; } - template - template - void MappingInfo::MappingInfoDependent::print_memory_consumption - (StreamType &out, - const SizeInfo &size_info) const + template + template + bool + FPArrayComparator:: + operator ()(const Tensor<1,dim,Tensor<1,VectorizedArray::n_array_elements,Number> > &t1, + const Tensor<1,dim,Tensor<1,VectorizedArray::n_array_elements,Number> > &t2) const { - // print_memory_statistics involves global communication, so we can - // disable the check here only if no processor has any such data -#ifdef DEAL_II_WITH_MPI - unsigned int general_size_glob = 0, general_size_loc = jacobians.size(); - int ierr = MPI_Allreduce (&general_size_loc, &general_size_glob, 1, - MPI_UNSIGNED, MPI_MAX, size_info.communicator); - AssertThrowMPI (ierr); -#else - unsigned int general_size_glob = jacobians.size(); -#endif - if (general_size_glob > 0) - { - out << " Memory Jacobian data: "; - size_info.print_memory_statistics - (out, MemoryConsumption::memory_consumption (jacobians) + - MemoryConsumption::memory_consumption (JxW_values)); - out << " Memory second derivative data: "; - size_info.print_memory_statistics - (out,MemoryConsumption::memory_consumption (jacobians_grad_diag) + - MemoryConsumption::memory_consumption (jacobians_grad_upper)); - } - -#ifdef DEAL_II_WITH_MPI - unsigned int quad_size_glob = 0, quad_size_loc = quadrature_points.size(); - ierr = MPI_Allreduce (&quad_size_loc, &quad_size_glob, 1, MPI_UNSIGNED, - MPI_MAX, size_info.communicator); - AssertThrowMPI (ierr); -#else - unsigned int quad_size_glob = quadrature_points.size(); -#endif - if (quad_size_glob > 0) - { - out << " Memory quadrature points: "; - size_info.print_memory_statistics - (out, MemoryConsumption::memory_consumption (rowstart_q_points) + - MemoryConsumption::memory_consumption (quadrature_points)); - } + for (unsigned int d=0; d::n_array_elements; ++k) + if (t1[d][k] < t2[d][k] - tolerance) + return true; + else if (t1[d][k] > t2[d][k] + tolerance) + return false; + return false; } - template - template - void MappingInfo::print_memory_consumption(StreamType &out, - const SizeInfo &size_info) const + template + template + bool + FPArrayComparator:: + operator ()(const Tensor<2,dim,Tensor<1,VectorizedArray::n_array_elements,Number> > &t1, + const Tensor<2,dim,Tensor<1,VectorizedArray::n_array_elements,Number> > &t2) const { - out << " Cell types: "; - size_info.print_memory_statistics - (out, MemoryConsumption::memory_consumption (cell_type)); - out << " Memory transformations compr: "; - size_info.print_memory_statistics - (out, MemoryConsumption::memory_consumption (affine_data) + - MemoryConsumption::memory_consumption (cartesian_data)); - for (unsigned int j=0; j::n_array_elements; ++k) + if (t1[d][e][k] < t2[d][e][k] - tolerance) + return true; + else if (t1[d][e][k] > t2[d][e][k] + tolerance) + return false; + return false; } } // end of namespace MatrixFreeFunctions diff --git a/include/deal.II/matrix_free/matrix_free.h b/include/deal.II/matrix_free/matrix_free.h index f94a1bf865..4d3daa83ed 100644 --- a/include/deal.II/matrix_free/matrix_free.h +++ b/include/deal.II/matrix_free/matrix_free.h @@ -1397,9 +1397,8 @@ unsigned int MatrixFree::get_n_q_points(const unsigned int quad_index, const unsigned int active_fe_index) const { - AssertIndexRange (quad_index, - mapping_info.mapping_data_gen.size()); - return mapping_info.mapping_data_gen[quad_index].n_q_points[active_fe_index]; + AssertIndexRange (quad_index, mapping_info.cell_data.size()); + return mapping_info.cell_data[quad_index].descriptor[active_fe_index].n_q_points; } @@ -1422,9 +1421,8 @@ unsigned int MatrixFree::get_n_q_points_face(const unsigned int quad_index, const unsigned int active_fe_index) const { - AssertIndexRange (quad_index, - mapping_info.mapping_data_gen.size()); - return mapping_info.mapping_data_gen[quad_index].n_q_points_face[active_fe_index]; + AssertIndexRange (quad_index, mapping_info.face_data.size()); + return mapping_info.face_data[quad_index].descriptor[active_fe_index].n_q_points; } @@ -1475,9 +1473,8 @@ const Quadrature & MatrixFree::get_quadrature (const unsigned int quad_index, const unsigned int active_fe_index) const { - AssertIndexRange (quad_index, mapping_info.mapping_data_gen.size()); - return mapping_info.mapping_data_gen[quad_index]. - quadrature[active_fe_index]; + AssertIndexRange (quad_index, mapping_info.cell_data.size()); + return mapping_info.cell_data[quad_index].descriptor[active_fe_index].quadrature; } @@ -1488,9 +1485,8 @@ const Quadrature & MatrixFree::get_face_quadrature (const unsigned int quad_index, const unsigned int active_fe_index) const { - AssertIndexRange (quad_index, mapping_info.mapping_data_gen.size()); - return mapping_info.mapping_data_gen[quad_index]. - face_quadrature[active_fe_index]; + AssertIndexRange (quad_index, mapping_info.face_data.size()); + return mapping_info.face_data[quad_index].descriptor[active_fe_index].quadrature; } diff --git a/include/deal.II/matrix_free/matrix_free.templates.h b/include/deal.II/matrix_free/matrix_free.templates.h index 22b8eade93..0864dc8f02 100644 --- a/include/deal.II/matrix_free/matrix_free.templates.h +++ b/include/deal.II/matrix_free/matrix_free.templates.h @@ -1034,73 +1034,6 @@ namespace internal boundary_cells_start = boundary_cells_end = n_macro_cells; } - - - /* ------------------------------------------------------------------ */ - - template - FPArrayComparator::FPArrayComparator (const Number scaling) - : - tolerance (scaling * std::numeric_limits::epsilon() * 1024.) - {} - - - - template - bool - FPArrayComparator::operator() (const std::vector &v1, - const std::vector &v2) const - { - const unsigned int s1 = v1.size(), s2 = v2.size(); - if (s1 < s2) - return true; - else if (s1 > s2) - return false; - else - for (unsigned int i=0; i v2[i] + tolerance) - return false; - return false; - } - - - - template - template - bool - FPArrayComparator:: - operator ()(const Tensor<1,dim,Tensor<1,VectorizedArray::n_array_elements,Number> > &t1, - const Tensor<1,dim,Tensor<1,VectorizedArray::n_array_elements,Number> > &t2) const - { - for (unsigned int d=0; d::n_array_elements; ++k) - if ((t1)[d][k] < (t2)[d][k] - tolerance) - return true; - else if ((t1)[d][k] > (t2)[d][k] + tolerance) - return false; - return false; - } - - - - template - template - bool - FPArrayComparator:: - operator ()(const Tensor<2,dim,Tensor<1,VectorizedArray::n_array_elements,Number> > &t1, - const Tensor<2,dim,Tensor<1,VectorizedArray::n_array_elements,Number> > &t2) const - { - for (unsigned int d=0; d::n_array_elements; ++k) - if ((t1)[d][e][k] < (t2)[d][e][k] - tolerance) - return true; - else if ((t1)[d][e][k] > (t2)[d][e][k] + tolerance) - return false; - return false; - } } } diff --git a/source/matrix_free/matrix_free.inst.in b/source/matrix_free/matrix_free.inst.in index a79ed61004..b15a166af6 100644 --- a/source/matrix_free/matrix_free.inst.in +++ b/source/matrix_free/matrix_free.inst.in @@ -32,6 +32,11 @@ for (deal_II_dimension : DIMENSIONS) template struct internal::MatrixFreeFunctions::MappingInfo; template struct internal::MatrixFreeFunctions::MappingInfo; + template struct internal::MatrixFreeFunctions::MappingInfoStorage; + template struct internal::MatrixFreeFunctions::MappingInfoStorage; + template struct internal::MatrixFreeFunctions::MappingInfoStorage; + template struct internal::MatrixFreeFunctions::MappingInfoStorage; + #ifndef DEAL_II_MSVC template void diff --git a/tests/matrix_free/compress_mapping.cc b/tests/matrix_free/compress_mapping.cc index 777fa0f01e..1507095666 100644 --- a/tests/matrix_free/compress_mapping.cc +++ b/tests/matrix_free/compress_mapping.cc @@ -79,14 +79,15 @@ void test () mf.reinit (dof, constraints, quad, data); const unsigned int n_macro_cells = mf.n_macro_cells(); - const unsigned int n_cartesian = mf.get_mapping_info().cartesian_data.size(); - const unsigned int n_affine = mf.get_mapping_info().affine_data.size(); - const unsigned int n_general = mf.get_mapping_info().mapping_data_gen[0].rowstart_jacobians.size()-1; + std::vector n_cell_types(4, 0); + for (unsigned int i=0; i 0, ExcInternalError()); + Assert(mf.get_mapping_info().cell_data[0].jacobians[0].size() < + (n_cell_types[3]*quad.size()+n_macro_cells-n_cell_types[3]), + ExcInternalError()); deallog << "OK" << std::endl; } @@ -115,15 +116,14 @@ void test_cube () mf.reinit (dof, constraints, quad, data); const unsigned int n_macro_cells = mf.n_macro_cells(); - const unsigned int n_cartesian = mf.get_mapping_info().cartesian_data.size(); - const unsigned int n_affine = mf.get_mapping_info().affine_data.size(); - const unsigned int n_general = mf.get_mapping_info().mapping_data_gen[0].rowstart_jacobians.size()-1; + std::vector n_cell_types(4, 0); + for (unsigned int i=0; i 1, ExcInternalError()); deallog << "OK" << std::endl; } @@ -160,15 +160,14 @@ void test_parallelogram () mf.reinit (dof, constraints, quad, data); const unsigned int n_macro_cells = mf.n_macro_cells(); - const unsigned int n_cartesian = mf.get_mapping_info().cartesian_data.size(); - const unsigned int n_affine = mf.get_mapping_info().affine_data.size(); - const unsigned int n_general = mf.get_mapping_info().mapping_data_gen[0].rowstart_jacobians.size()-1; + std::vector n_cell_types(4, 0); + for (unsigned int i=0; i 1, ExcInternalError()); deallog << "OK" << std::endl; }