From 9fb6bd59940f4cc3e81908f06a14bb57b30729ce Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Mon, 23 Nov 2015 13:31:20 -0600 Subject: [PATCH] Do some work already in get_data(). Do not delay until fill_fe_values(), if something can already be done early on. --- include/deal.II/fe/fe_poly.h | 110 +++++++++++++++---------- include/deal.II/fe/fe_poly.templates.h | 14 ++-- source/fe/fe_poly.cc | 29 +++---- 3 files changed, 88 insertions(+), 65 deletions(-) diff --git a/include/deal.II/fe/fe_poly.h b/include/deal.II/fe/fe_poly.h index 134b5ef3a0..01db90689f 100644 --- a/include/deal.II/fe/fe_poly.h +++ b/include/deal.II/fe/fe_poly.h @@ -223,7 +223,7 @@ protected: get_data(const UpdateFlags update_flags, const Mapping &/*mapping*/, const Quadrature &quadrature, - dealii::internal::FEValues::FiniteElementRelatedData &/*output_data*/) const + dealii::internal::FEValues::FiniteElementRelatedData &output_data) const { // generate a new data object and // initialize some fields @@ -232,51 +232,58 @@ protected: const unsigned int n_q_points = quadrature.size(); - // some scratch arrays - std::vector values(0); - std::vector > grads(0); - std::vector > grad_grads(0); - std::vector > third_derivatives(0); - std::vector > fourth_derivatives(0); - - // initialize fields only if really - // necessary. otherwise, don't - // allocate memory - if (update_flags & update_values) - { - values.resize (this->dofs_per_cell); - data->shape_values.resize (this->dofs_per_cell, - std::vector (n_q_points)); - } + // initialize some scratch arrays. we need them for the underlying + // polynomial to put the values and derivatives of shape functions + // to put there, depending on what the user requested + std::vector values(update_flags & update_values ? + this->dofs_per_cell : 0); + std::vector > grads(update_flags & update_gradients ? + this->dofs_per_cell : 0); + std::vector > grad_grads(update_flags & update_hessians ? + this->dofs_per_cell : 0); + std::vector > third_derivatives(update_flags & update_3rd_derivatives ? + this->dofs_per_cell : 0); + std::vector > fourth_derivatives; // won't be needed, so leave empty + + // now also initialize fields the fields of this class's own + // temporary storage, depending on what we need for the given + // update flags. + // + // there is one exception from the rule: if we are dealing with + // cells (i.e., if this function is not called via + // get_(sub)face_data()), then we can already store things in the + // final location where FEValues::reinit() later wants to see + // things. we then don't need the intermediate space. we determine + // whether we are on a cell by asking whether the number of + // elements in the output array equals the number of quadrature + // points (yes, it's a cell) or not (because in that case the + // number of quadrature points we use here equals the number of + // quadrature points summed over *all* faces or subfaces, whereas + // the number of output slots equals the number of quadrature + // points on only *one* face) + if ((update_flags & update_values) + && + !((output_data.shape_values.n_rows() > 0) + && + (output_data.shape_values.n_cols() == n_q_points))) + data->shape_values.resize (this->dofs_per_cell, + std::vector (n_q_points)); if (update_flags & update_gradients) - { - grads.resize (this->dofs_per_cell); - data->shape_gradients.resize (this->dofs_per_cell, - std::vector > (n_q_points)); - } + data->shape_gradients.resize (this->dofs_per_cell, + std::vector > (n_q_points)); if (update_flags & update_hessians) - { - grad_grads.resize (this->dofs_per_cell); - data->shape_hessians.resize (this->dofs_per_cell, - std::vector > (n_q_points)); - } + data->shape_hessians.resize (this->dofs_per_cell, + std::vector > (n_q_points)); if (update_flags & update_3rd_derivatives) - { - third_derivatives.resize (this->dofs_per_cell); - data->shape_3rd_derivatives.resize (this->dofs_per_cell, - std::vector > (n_q_points)); - } - - // next already fill those fields - // of which we have information by - // now. note that the shape - // gradients are only those on the - // unit cell, and need to be - // transformed when visiting an - // actual cell + data->shape_3rd_derivatives.resize (this->dofs_per_cell, + std::vector > (n_q_points)); + + // next already fill those fields of which we have information by + // now. note that the shape gradients are only those on the unit + // cell, and need to be transformed when visiting an actual cell if (update_flags & (update_values | update_gradients | update_hessians | update_3rd_derivatives) ) for (unsigned int i=0; idofs_per_cell; ++k) - data->shape_values[k][i] = values[k]; - + if (output_data.shape_values.n_rows() > 0) + { + if (output_data.shape_values.n_cols() == n_q_points) + for (unsigned int k=0; kdofs_per_cell; ++k) + output_data.shape_values[k][i] = values[k]; + else + for (unsigned int k=0; kdofs_per_cell; ++k) + data->shape_values[k][i] = values[k]; + } + + // for everything else, derivatives need to be transformed, + // so we write them into our scratch space and only later + // copy stuff into where FEValues wants it if (update_flags & update_gradients) for (unsigned int k=0; kdofs_per_cell; ++k) data->shape_gradients[k][i] = grads[k]; diff --git a/include/deal.II/fe/fe_poly.templates.h b/include/deal.II/fe/fe_poly.templates.h index c55bf59053..bc2bf1d9d6 100644 --- a/include/deal.II/fe/fe_poly.templates.h +++ b/include/deal.II/fe/fe_poly.templates.h @@ -254,11 +254,9 @@ fill_fe_values (const typename Triangulation::cell_iterator &, const UpdateFlags flags(fe_data.update_each); - if (flags & update_values) - for (unsigned int k=0; kdofs_per_cell; ++k) - for (unsigned int i=0; idofs_per_cell; ++k) mapping.transform (fe_data.shape_gradients[k], @@ -329,6 +327,9 @@ fill_fe_face_values (const typename Triangulation::cell_iterator const UpdateFlags flags(fe_data.update_each); + // transform gradients and higher derivatives. we also have to copy + // the values (unlike in the case of fill_fe_values()) since + // we need to take into account the offsets if (flags & update_values) for (unsigned int k=0; kdofs_per_cell; ++k) for (unsigned int i=0; i::cell_iterato const UpdateFlags flags(fe_data.update_each); + // transform gradients and higher derivatives. we also have to copy + // the values (unlike in the case of fill_fe_values()) since + // we need to take into account the offsets if (flags & update_values) for (unsigned int k=0; kdofs_per_cell; ++k) for (unsigned int i=0; i::cell_iterator &, for (unsigned int k=0; kdofs_per_cell; ++k) { - if (fe_data.update_each & update_values) - for (unsigned int i=0; i::cell_iterator &, for (unsigned int k=0; kdofs_per_cell; ++k) { - if (fe_data.update_each & update_values) - for (unsigned int i=0; i::cell_iterator &, for (unsigned int k=0; kdofs_per_cell; ++k) { - if (fe_data.update_each & update_values) - for (unsigned int i=0; i::cell_iterator &, for (unsigned int k=0; kdofs_per_cell; ++k) { - if (fe_data.update_each & update_values) - for (unsigned int i=0; i