From 59fcaf9b7338c3c6363937775b56e1ef8a3a1ef1 Mon Sep 17 00:00:00 2001 From: Martin Kronbichler Date: Sat, 28 Oct 2017 17:04:44 +0200 Subject: [PATCH] Fix dofs_per_cell for FEEvaluation. --- include/deal.II/matrix_free/fe_evaluation.h | 98 +++++++++++++-------- include/deal.II/matrix_free/operators.h | 11 +-- tests/matrix_free/assemble_matrix_02.cc | 2 +- 3 files changed, 66 insertions(+), 45 deletions(-) diff --git a/include/deal.II/matrix_free/fe_evaluation.h b/include/deal.II/matrix_free/fe_evaluation.h index 45aeb0762a..f1121778ef 100644 --- a/include/deal.II/matrix_free/fe_evaluation.h +++ b/include/deal.II/matrix_free/fe_evaluation.h @@ -304,10 +304,13 @@ public: /** * Return the value stored for the local degree of freedom with index @p * dof. If the object is vector-valued, a vector-valued return argument is - * given. Note that when vectorization is enabled, values from several cells - * are grouped together. If @p set_dof_values was called last, the value - * corresponds to the one set there. If @p integrate was called last, it - * instead corresponds to the value of the integrated function with the test + * given. Thus, the argument @p dof can at most run until @p + * dofs_per_component rather than @p dofs_per_cell since the different + * components of a vector-valued FE are return together. Note that when + * vectorization is enabled, values from several cells are grouped + * together. If @p set_dof_values was called last, the value corresponds to + * the one set there. If @p integrate was called last, it instead + * corresponds to the value of the integrated function with the test * function of the given index. * * Note that the derived class FEEvaluationAccess overloads this operation @@ -594,7 +597,7 @@ public: * Returns an ArrayView to internal memory for temporary use. Note that some * of this memory is overwritten during evaluate() and integrate() calls so * do not assume it to be stable over those calls. The maximum size you can - * write into is (n_components+2)*dofs_per_cell+2*n_q_points. + * write into is 3*dofs_per_cell+2*n_q_points. */ ArrayView > get_scratch_data() const; @@ -1910,7 +1913,9 @@ public: static const unsigned int dimension = dim; static const unsigned int n_components = n_components_; static const unsigned int static_n_q_points = Utilities::fixed_int_power::value; + static const unsigned int static_dofs_per_component = Utilities::fixed_int_power::value; static const unsigned int tensor_dofs_per_cell = Utilities::fixed_int_power::value; + static const unsigned int static_dofs_per_cell = n_components * Utilities::fixed_int_power::value; /** * Constructor. Takes all data stored in MatrixFree. If applied to problems @@ -2024,16 +2029,26 @@ public: quadrature_point (const unsigned int q_point) const; /** - * The number of scalar degrees of freedom on the cell. Usually close to - * tensor_dofs_per_cell, but depends on the actual element selected and is - * thus not static. + * The number of degrees of freedom of a single component on the cell for + * the underlying evaluation object. Usually close to + * static_dofs_per_component, but the number depends on the actual element + * selected and is thus not static. + */ + const unsigned int dofs_per_component; + + /** + * The number of degrees of freedom on the cell accumulated over all + * components in the current evaluation object. Usually close to + * static_dofs_per_component*n_components=static_dofs_per_cell, but the + * number depends on the actual element selected and is thus not static. */ const unsigned int dofs_per_cell; /** - * The number of scalar degrees of freedom on the cell. Usually close to - * tensor_dofs_per_cell, but depends on the actual element selected and is - * thus not static. + * The number of quadrature points in use for the current evaluation + * object. It equals static_n_q_points in case the element degree is not set + * to -1, i.e., the worker kernels can be set according to the templates @p + * fe_degree and @p n_q_points_1d rather than at run time. */ const unsigned int n_q_points; @@ -2055,16 +2070,16 @@ namespace internal // a helper function to compute the number of DoFs of a DGP element at compile // time, depending on the degree template - struct DGP_dofs_per_cell + struct DGP_dofs_per_component { // this division is always without remainder static const unsigned int value = - (DGP_dofs_per_cell::value * (degree+dim)) / dim; + (DGP_dofs_per_component::value * (degree+dim)) / dim; }; // base specialization: 1d elements have 'degree+1' degrees of freedom template - struct DGP_dofs_per_cell<1,degree> + struct DGP_dofs_per_component<1,degree> { static const unsigned int value = degree+1; }; @@ -2132,8 +2147,8 @@ FEEvaluationBase ExcNotInitialized()); AssertDimension (matrix_info->get_size_info().vectorization_length, VectorizedArray::n_array_elements); - AssertDimension (data->dofs_per_cell, - dof_info->dofs_per_cell[active_fe_index]/n_fe_components); + AssertDimension (data->dofs_per_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]); Assert (n_fe_components == 1 || @@ -2377,33 +2392,33 @@ FEEvaluationBase { Assert(scratch_data_array != nullptr, ExcInternalError()); - const unsigned int tensor_dofs_per_cell = + const unsigned int tensor_dofs_per_component = Utilities::fixed_power(this->data->fe_degree+1); - const unsigned int dofs_per_cell = this->data->dofs_per_cell; + const unsigned int dofs_per_component = this->data->dofs_per_cell; const unsigned int n_quadrature_points = this->data->n_q_points; - const unsigned int shift = std::max(tensor_dofs_per_cell+1, dofs_per_cell)* - (n_components_+2) + 2 * n_quadrature_points; - const unsigned int allocated_size = shift + n_components_ * dofs_per_cell + const unsigned int shift = std::max(tensor_dofs_per_component+1, dofs_per_component)* + n_components_*3 + 2 * n_quadrature_points; + const unsigned int allocated_size = shift + n_components_ * dofs_per_component + (n_components_*(dim*dim+2*dim+1)*n_quadrature_points); scratch_data_array->resize_fast(allocated_size); // set the pointers to the correct position in the data array for (unsigned int c=0; cvalues_dofs[c] = scratch_data_array->begin() + c*dofs_per_cell; + this->values_dofs[c] = scratch_data_array->begin() + c*dofs_per_component; this->values_quad[c] = scratch_data_array->begin() + - n_components*dofs_per_cell+c*n_quadrature_points; + n_components*dofs_per_component+c*n_quadrature_points; for (unsigned int d=0; dgradients_quad[c][d] = scratch_data_array->begin() + - n_components*(dofs_per_cell+n_quadrature_points) + + n_components*(dofs_per_component+n_quadrature_points) + (c*dim+d)*n_quadrature_points; for (unsigned int d=0; d<(dim*dim+dim)/2; ++d) this->hessians_quad[c][d] = scratch_data_array->begin() + - n_components*((dim+1)*n_quadrature_points + dofs_per_cell) + + n_components*((dim+1)*n_quadrature_points + dofs_per_component) + (c*(dim*dim+dim)+d)*n_quadrature_points; } - scratch_data = scratch_data_array->begin() + n_components_ * dofs_per_cell + scratch_data = scratch_data_array->begin() + n_components_ * dofs_per_component + (n_components_*(dim*dim+2*dim+1)*n_quadrature_points); } @@ -3021,7 +3036,7 @@ FEEvaluationBase const std::pair *indicators_end = dof_info->end_indicators(cell); unsigned int ind_local = 0; - const unsigned int dofs_per_cell = this->data->dofs_per_cell; + const unsigned int dofs_per_component = this->data->dofs_per_cell; const unsigned int n_irreg_components_filled = dof_info->row_starts[cell][2]; const bool at_irregular_cell = n_irreg_components_filled > 0; @@ -3041,7 +3056,7 @@ FEEvaluationBase "read_dof_values and distribute_local_to_global.")); const unsigned int n_local_dofs = - VectorizedArray::n_array_elements * dofs_per_cell; + VectorizedArray::n_array_elements * dofs_per_component; for (unsigned int comp=0; comp // vectorization loop AssertDimension (dof_info->end_indices(cell)-dof_indices, static_cast(n_local_dofs)); - for (unsigned int j=0, ind=0; j::n_array_elements) + for (unsigned int j=0, ind=0; j::n_array_elements) for (unsigned int comp=0; comp internal::check_vector_compatibility (*src[0], *dof_info); Assert (n_fe_components == n_components_, ExcNotImplemented()); const unsigned int n_local_dofs = - dofs_per_cell*VectorizedArray::n_array_elements * n_components; + dofs_per_component*VectorizedArray::n_array_elements * n_components; Number *local_data = const_cast(&values_dofs[0][0][0]); if (at_irregular_cell == false) @@ -3248,7 +3263,7 @@ FEEvaluationBase AssertDimension (dof_info->end_indices(cell)-dof_indices, static_cast(n_local_dofs)); for (unsigned int comp=0, ind=0; comp::n_array_elements) + for (unsigned int j=0; j::n_array_elements) operation.process_dof_gather(dof_indices+ind, *src[0], values_dofs[comp][j], std::integral_constant::value>()); @@ -3447,7 +3462,7 @@ FEEvaluationBase // iterates over the elements of index_local_to_global and dof_indices // points to the global indices stored in index_local_to_global const unsigned int *dof_indices = dof_info->begin_indices_plain(cell); - const unsigned int dofs_per_cell = this->data->dofs_per_cell; + const unsigned int dofs_per_component = this->data->dofs_per_cell; const unsigned int n_irreg_components_filled = dof_info->row_starts[cell][2]; const bool at_irregular_cell = n_irreg_components_filled > 0; @@ -3467,7 +3482,7 @@ FEEvaluationBase "read_dof_values_plain.")); const unsigned int n_local_dofs = - VectorizedArray::n_array_elements * dofs_per_cell; + VectorizedArray::n_array_elements * dofs_per_component; for (unsigned int comp=0; comp internal::check_vector_compatibility (*src[0], *dof_info); Assert (n_fe_components == n_components_, ExcNotImplemented()); const unsigned int n_local_dofs = - dofs_per_cell * VectorizedArray::n_array_elements * n_components; + dofs_per_component * VectorizedArray::n_array_elements * n_components; Number *local_src_number = &values_dofs[0][0][0]; if (at_irregular_cell == false) { @@ -5084,7 +5099,8 @@ FEEvaluation const unsigned int quad_no) : BaseClass (data_in, fe_no, quad_no, fe_degree, static_n_q_points), - dofs_per_cell (this->data->dofs_per_cell), + dofs_per_component (this->data->dofs_per_cell), + dofs_per_cell (this->data->dofs_per_cell*n_components_), n_q_points (this->data->n_q_points) { check_template_arguments(fe_no, 0); @@ -5105,7 +5121,8 @@ FEEvaluation BaseClass (mapping, fe, quadrature, update_flags, first_selected_component, static_cast*>(nullptr)), - dofs_per_cell (this->data->dofs_per_cell), + dofs_per_component (this->data->dofs_per_cell), + dofs_per_cell (this->data->dofs_per_cell*n_components_), n_q_points (this->data->n_q_points) { check_template_arguments(numbers::invalid_unsigned_int, 0); @@ -5125,7 +5142,8 @@ FEEvaluation BaseClass (StaticMappingQ1::mapping, fe, quadrature, update_flags, first_selected_component, static_cast*>(nullptr)), - dofs_per_cell (this->data->dofs_per_cell), + dofs_per_component (this->data->dofs_per_cell), + dofs_per_cell (this->data->dofs_per_cell*n_components_), n_q_points (this->data->n_q_points) { check_template_arguments(numbers::invalid_unsigned_int, 0); @@ -5146,7 +5164,8 @@ FEEvaluation fe, other.mapped_geometry->get_quadrature(), other.mapped_geometry->get_fe_values().get_update_flags(), first_selected_component, &other), - dofs_per_cell (this->data->dofs_per_cell), + dofs_per_component (this->data->dofs_per_cell), + dofs_per_cell (this->data->dofs_per_cell*n_components_), n_q_points (this->data->n_q_points) { check_template_arguments(numbers::invalid_unsigned_int, 0); @@ -5161,7 +5180,8 @@ FEEvaluation ::FEEvaluation (const FEEvaluation &other) : BaseClass (other), - dofs_per_cell (this->data->dofs_per_cell), + dofs_per_component (this->data->dofs_per_cell), + dofs_per_cell (this->data->dofs_per_cell*n_components_), n_q_points (this->data->n_q_points) { check_template_arguments(numbers::invalid_unsigned_int, 0); diff --git a/include/deal.II/matrix_free/operators.h b/include/deal.II/matrix_free/operators.h index cc65ed759c..2c20b1f6c9 100644 --- a/include/deal.II/matrix_free/operators.h +++ b/include/deal.II/matrix_free/operators.h @@ -1781,17 +1781,18 @@ namespace MatrixFreeOperators for (unsigned int cell=cell_range.first; cell local_diagonal_vector[phi.tensor_dofs_per_cell]; - for (unsigned int i=0; i local_diagonal_vector[phi.static_dofs_per_cell]; + for (unsigned int i=0; i(); phi.begin_dof_values()[i] = 1.; do_operation_on_cell(phi,cell); local_diagonal_vector[i] = phi.begin_dof_values()[i]; } - for (unsigned int i=0; i &dof) phi_u.reinit(cell); phi_p.reinit(cell); - const unsigned int dofs_per_cell_u = phi_u.dofs_per_cell*dim; + const unsigned int dofs_per_cell_u = phi_u.dofs_per_cell; const unsigned int dofs_per_cell_p = phi_p.dofs_per_cell; for (unsigned int i=0; i::n_array_elements) -- 2.39.5