From f89865812ce7166850c215d71e01f2173a4947f5 Mon Sep 17 00:00:00 2001 From: bangerth Date: Thu, 11 Dec 2008 19:01:53 +0000 Subject: [PATCH] Modify the FEValuesViews classes so that they compute a bunch of things already at construction time. This should make them much faster when you use them to evaluate stuff. On the other hand, it makes construction much more expensive and so the previous approach of creating them on the fly every time someone wrote fe_values[velocities] doesn't work any more. Rather, create all possible views objects at construction time of the FEValues object and simply return a reference to one of them when a view is requested. git-svn-id: https://svn.dealii.org/trunk@17913 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/deal.II/include/fe/fe_values.h | 752 ++++++++----------------- deal.II/deal.II/source/fe/fe_values.cc | 245 +++++++- 2 files changed, 457 insertions(+), 540 deletions(-) diff --git a/deal.II/deal.II/include/fe/fe_values.h b/deal.II/deal.II/include/fe/fe_values.h index 5798daaef2..28d39c7b53 100644 --- a/deal.II/deal.II/include/fe/fe_values.h +++ b/deal.II/deal.II/include/fe/fe_values.h @@ -184,6 +184,12 @@ namespace FEValuesViews */ typedef Tensor<2,spacedim> hessian_type; + /** + * Default constructor. Creates an + * invalid object. + */ + Scalar (); + /** * Constructor for an object that * represents a single scalar component @@ -194,6 +200,14 @@ namespace FEValuesViews Scalar (const FEValuesBase &fe_values_base, const unsigned int component); + /** + * Copy operator. This is not a + * lightweight object so we don't allow + * copying and generate an exception if + * this function is called. + */ + Scalar & operator= (const Scalar &); + /** * Return the value of the vector * component selected by this view, for @@ -241,6 +255,42 @@ namespace FEValuesViews * object. */ const unsigned int component; + + /** + * For each shape function, store + * whether it is primitive. + */ + std::vector is_primitive; + + /** + * For each shape function, store + * whether the selected vector + * component may be nonzero. For + * primitive shape functions we know + * for sure whether a certain scalar + * component of a given shape function + * is nonzero, whereas for + * non-primitive shape functions this + * may not be entirely clear (e.g. for + * RT elements it depends on the shape + * of a cell). + */ + Table<1,bool> is_nonzero_shape_function_component; + + /** + * For each shape function, store the + * row index within the shape_values, + * shape_gradients, and shape_hessians + * tables (the column index is the + * quadrature point index). If the + * shape function is primitive, then we + * can get this information from the + * shape_function_to_row_table of the + * FEValues object; otherwise, we have + * to work a bit harder to compute this + * information. + */ + Table<1,unsigned int> row_index; }; @@ -308,6 +358,12 @@ namespace FEValuesViews */ typedef Tensor<3,dim> hessian_type; + /** + * Default constructor. Creates an + * invalid object. + */ + Vector (); + /** * Constructor for an object that * represents dim components of a @@ -322,7 +378,15 @@ namespace FEValuesViews */ Vector (const FEValuesBase &fe_values_base, const unsigned int first_vector_component); - + + /** + * Copy operator. This is not a + * lightweight object so we don't allow + * copying and generate an exception if + * this function is called. + */ + Vector & operator= (const Vector &); + /** * Return the value of the vector * components selected by this view, @@ -407,10 +471,77 @@ namespace FEValuesViews * FEValuesBase object. */ const unsigned int first_vector_component; + + /** + * For each pair (shape + * function,component within vector), + * store whether the selected vector + * component may be nonzero. For + * primitive shape functions we know + * for sure whether a certain scalar + * component of a given shape function + * is nonzero, whereas for + * non-primitive shape functions this + * may not be entirely clear (e.g. for + * RT elements it depends on the shape + * of a cell). + */ + Table<2,bool> is_nonzero_shape_function_component; + + /** + * For each pair (shape function, + * component within vector), store the + * row index within the shape_values, + * shape_gradients, and shape_hessians + * tables (the column index is the + * quadrature point index). If the + * shape function is primitive, then we + * can get this information from the + * shape_function_to_row_table of the + * FEValues object; otherwise, we have + * to work a bit harder to compute this + * information. + */ + Table<2,unsigned int> row_index; }; } +namespace internal +{ + namespace FEValuesViews + { + /** + * A class objects of which store a + * collection of FEValuesViews::Scalar, + * FEValuesViews::Vector, etc object. The + * FEValuesBase class uses it to generate + * all possible Views classes upon + * construction time; we do this at + * construction time since the Views + * classes cache some information and are + * therefore relatively expensive to + * create. + */ + template + struct Cache + { + /** + * Caches for scalar and + * vector-valued views. + */ + std::vector > scalars; + std::vector > vectors; + + /** + * Constructor. + */ + Cache (const FEValuesBase &fe_values); + }; + } +} + + /*!@addtogroup feaccess */ @@ -866,7 +997,7 @@ template * FEValuesViews and in particular * in the @ref vector_valued module. */ - FEValuesViews::Scalar + const FEValuesViews::Scalar & operator[] (const FEValuesExtractors::Scalar &scalar) const; /** @@ -879,7 +1010,7 @@ template * the namespace FEValuesViews and in particular * in the @ref vector_valued module. */ - FEValuesViews::Vector + const FEValuesViews::Vector & operator[] (const FEValuesExtractors::Vector &vector) const; //@} @@ -2202,7 +2333,7 @@ template */ const CI cell; }; - + /** * Implementation of a derived @@ -2511,6 +2642,12 @@ template */ FEValuesBase & operator= (const FEValuesBase &); + /** + * A cache for all possible FEValuesViews + * objects. + */ + internal::FEValuesViews::Cache fe_values_views_cache; + /** * Make the view classes friends of this * class, since they access internal @@ -2532,7 +2669,7 @@ template * @author Wolfgang Bangerth, 1998, Guido Kanschat, 2001 */ template - class FEValues : public FEValuesBase +class FEValues : public FEValuesBase { public: /** @@ -3252,20 +3389,6 @@ namespace FEValuesExtractors namespace FEValuesViews { - template - inline - Scalar::Scalar (const FEValuesBase &fe_values, - const unsigned int component) - : - fe_values (fe_values), - component (component) - { - Assert (component < fe_values.fe->n_components(), - ExcIndexRange(component, 0, fe_values.fe->n_components())); - } - - - template inline typename Scalar::value_type @@ -3281,37 +3404,15 @@ namespace FEValuesViews // an adaptation of the // FEValuesBase::shape_value_component // function except that here we know the - // component as fixed. see the - // comments there - // - // we can do away with some of the - // assertions since they are already - // taken care of in the constructor of - // this class - if (fe_values.fe->is_primitive() || - fe_values.fe->is_primitive(shape_function)) - { - if (component == - fe_values.fe->system_to_component_index(shape_function).first) - return fe_values.shape_values(fe_values.shape_function_to_row_table[shape_function],q_point); - else - return 0; - } + // component as fixed and we have + // pre-computed and cached a bunch of + // information. see the comments there + if (is_nonzero_shape_function_component[shape_function]) + return fe_values.shape_values(row_index[shape_function],q_point); else - { - if (fe_values.fe->get_nonzero_components(shape_function)[component] == false) - return 0; - - const unsigned int - row = (fe_values.shape_function_to_row_table[shape_function] - + - std::count (fe_values.fe->get_nonzero_components(shape_function).begin(), - fe_values.fe->get_nonzero_components(shape_function).begin()+ - component, - true)); - return fe_values.shape_values(row, q_point); - } + return 0; } + @@ -3330,36 +3431,13 @@ namespace FEValuesViews // an adaptation of the // FEValuesBase::shape_grad_component // function except that here we know the - // component as fixed. see the - // comments there - // - // we can do away with the assertions - // since they are already taken care of - // in the constructor of this class - if (fe_values.fe->is_primitive() || - fe_values.fe->is_primitive(shape_function)) - { - if (component == - fe_values.fe->system_to_component_index(shape_function).first) - return fe_values.shape_gradients[fe_values. - shape_function_to_row_table[shape_function]][q_point]; - else - return gradient_type(); - } + // component as fixed and we have + // pre-computed and cached a bunch of + // information. see the comments there + if (is_nonzero_shape_function_component[shape_function]) + return fe_values.shape_gradients[row_index[shape_function]][q_point]; else - { - if (fe_values.fe->get_nonzero_components(shape_function)[component] == false) - return gradient_type(); - - const unsigned int - row = (fe_values.shape_function_to_row_table[shape_function] - + - std::count (fe_values.fe->get_nonzero_components(shape_function).begin(), - fe_values.fe->get_nonzero_components(shape_function).begin()+ - component, - true)); - return fe_values.shape_gradients[row][q_point]; - } + return gradient_type(); } @@ -3379,51 +3457,13 @@ namespace FEValuesViews // an adaptation of the // FEValuesBase::shape_grad_component // function except that here we know the - // component as fixed. see the - // comments there - // - // we can do away with the assertions - // since they are already taken care of - // in the constructor of this class - if (fe_values.fe->is_primitive() || - fe_values.fe->is_primitive(shape_function)) - { - if (component == - fe_values.fe->system_to_component_index(shape_function).first) - return fe_values.shape_hessians[fe_values. - shape_function_to_row_table[shape_function]][q_point]; - else - return hessian_type(); - } + // component as fixed and we have + // pre-computed and cached a bunch of + // information. see the comments there + if (is_nonzero_shape_function_component[shape_function]) + return fe_values.shape_hessians[row_index[shape_function]][q_point]; else - { - if (fe_values.fe->get_nonzero_components(shape_function)[component] == false) - return hessian_type(); - - const unsigned int - row = (fe_values.shape_function_to_row_table[shape_function] - + - std::count (fe_values.fe->get_nonzero_components(shape_function).begin(), - fe_values.fe->get_nonzero_components(shape_function).begin()+ - component, - true)); - return fe_values.shape_hessians[row][q_point]; - } - } - - - - template - inline - Vector::Vector (const FEValuesBase &fe_values, - const unsigned int first_vector_component) - : - fe_values (fe_values), - first_vector_component (first_vector_component) - { - Assert (first_vector_component+dim-1 < this->fe_values.fe->n_components(), - ExcIndexRange(first_vector_component+dim-1, 0, - this->fe_values.fe->n_components())); + return hessian_type(); } @@ -3440,53 +3480,19 @@ namespace FEValuesViews Assert (fe_values.update_flags & update_values, typename FVB::ExcAccessToUninitializedField()); - // compared to the scalar case above, we - // can save some work here because we - // know that we are querying a contiguous - // range of components + // same as for the scalar case except + // that we have one more index + // + // for primitive elements we could + // probably do even better than the loop + // below because we then know that only + // for one value of 'd' the + // 'if'-condition is true value_type return_value; - - if (fe_values.fe->is_primitive() || - fe_values.fe->is_primitive(shape_function)) - { - // if this is a primitive shape - // function then at most one element - // of the output vector is - // nonzero. find out if indeed one is - const unsigned int - nonzero_component - = fe_values.fe->system_to_component_index(shape_function).first; - - if ((nonzero_component >= first_vector_component) - && - (nonzero_component < first_vector_component + dim)) - return_value[nonzero_component - first_vector_component] - = fe_values.shape_values(fe_values. - shape_function_to_row_table[shape_function], - q_point); - } - else - { - unsigned int - row = (fe_values.shape_function_to_row_table[shape_function] - + - std::count (fe_values.fe->get_nonzero_components(shape_function).begin(), - fe_values.fe->get_nonzero_components(shape_function).begin() + - first_vector_component, - true)); - for (unsigned int d=0; dget_nonzero_components(shape_function)[first_vector_component+d] == - true) - { - return_value[d] = fe_values.shape_values(row, q_point); - - if ((d != dim-1) - && - (fe_values.fe->get_nonzero_components(shape_function)[first_vector_component+d] - == true)) - ++row; - } - } + for (unsigned int d=0; dis_primitive() || - fe_values.fe->is_primitive(shape_function)) - { - const unsigned int - nonzero_component - = fe_values.fe->system_to_component_index(shape_function).first; - - if ((nonzero_component >= first_vector_component) - && - (nonzero_component < first_vector_component + dim)) - return_value[nonzero_component - first_vector_component] - = fe_values.shape_gradients[fe_values. - shape_function_to_row_table[shape_function]][q_point]; - } - else - { - unsigned int - row = (fe_values.shape_function_to_row_table[shape_function] - + - std::count (fe_values.fe->get_nonzero_components(shape_function).begin(), - fe_values.fe->get_nonzero_components(shape_function).begin() + - first_vector_component, - true)); - for (unsigned int d=0; dget_nonzero_components(shape_function)[first_vector_component+d] == - true) - { - return_value[d] = fe_values.shape_gradients[row][q_point]; - - if ((d != dim-1) - && - (fe_values.fe->get_nonzero_components(shape_function)[first_vector_component+d] - == true)) - ++row; - } - } + for (unsigned int d=0; dis_primitive() || - fe_values.fe->is_primitive(shape_function)) - { - const unsigned int - nonzero_component - = fe_values.fe->system_to_component_index(shape_function).first; - - if ((nonzero_component >= first_vector_component) - && - (nonzero_component < first_vector_component + dim)) - return - fe_values.shape_gradients[fe_values. - shape_function_to_row_table[shape_function]][q_point] - [nonzero_component - first_vector_component]; - else - return 0; - } - else - { - unsigned int - row = (fe_values.shape_function_to_row_table[shape_function] - + - std::count (fe_values.fe->get_nonzero_components(shape_function).begin(), - fe_values.fe->get_nonzero_components(shape_function).begin() + - first_vector_component, - true)); - - double div = 0; - for (unsigned int d=0; dget_nonzero_components(shape_function)[first_vector_component+d] == - true) - { - div += fe_values.shape_gradients[row][q_point][d]; - - if ((d != dim-1) - && - (fe_values.fe->get_nonzero_components(shape_function)[first_vector_component+d] - == true)) - ++row; - } - return div; - } + // same as for the scalar case except + // that we have one more index + // + // for primitive elements we could + // probably do even better than the loop + // below because we then know that only + // for one value of 'd' the + // 'if'-condition is true + divergence_type return_value = 0; + for (unsigned int d=0; dis_primitive() || - fe_values.fe->is_primitive(shape_function)) - { - const unsigned int - nonzero_component - = fe_values.fe->system_to_component_index(shape_function).first; - - if ((nonzero_component >= first_vector_component) - && - (nonzero_component < first_vector_component + dim)) - return_value[nonzero_component - first_vector_component] - = fe_values.shape_hessians[fe_values. - shape_function_to_row_table[shape_function]][q_point]; - } - else - { - unsigned int - row = (fe_values.shape_function_to_row_table[shape_function] - + - std::count (fe_values.fe->get_nonzero_components(shape_function).begin(), - fe_values.fe->get_nonzero_components(shape_function).begin() + - first_vector_component, - true)); - for (unsigned int d=0; dget_nonzero_components(shape_function)[first_vector_component+d] == - true) - { - return_value[d] = fe_values.shape_hessians[row][q_point]; - - if ((d != dim-1) - && - (fe_values.fe->get_nonzero_components(shape_function)[first_vector_component+d] - == true)) - ++row; - } - } + for (unsigned int d=0; d - inline - Vector<1>::symmetric_gradient_type - Vector<1>::symmetric_gradient (const unsigned int shape_function, - const unsigned int q_point) const - { - // this function works like in the case - // above - Assert (shape_function < fe_values.fe->dofs_per_cell, - ExcIndexRange (shape_function, 0, fe_values.fe->dofs_per_cell)); - Assert (fe_values.update_flags & update_gradients, - FEValuesBase<1>::ExcAccessToUninitializedField()); - - if (fe_values.fe->is_primitive() || - fe_values.fe->is_primitive(shape_function)) - { - const unsigned int - nonzero_component - = fe_values.fe->system_to_component_index(shape_function).first; - - if (nonzero_component == first_vector_component) - { - // first get the one component of - // the nonsymmetrized gradient - // that is not zero - const Tensor<1,1> grad - = fe_values.shape_gradients[fe_values. - shape_function_to_row_table[shape_function]][q_point]; - - // then form a symmetric tensor - // out of it. note that access to - // individual elements of a - // SymmetricTensor object is - // fairly slow. if we implemented - // the following in the naive - // way, we would therefore incur - // a very significant penalty: a - // preliminary version of the - // Stokes tutorial program would - // slow down from 17 to 23 - // seconds because of this single - // function! - // - // consequently, we try to be a - // bit smarter by already laying - // out the data in the right - // format and creating a - // symmetric tensor of it - const double array[symmetric_gradient_type::n_independent_components] - = { grad[0] }; - return symmetric_gradient_type(array); - } - else - return symmetric_gradient_type(); - } - else - { - return symmetrize (gradient(shape_function, q_point)); - } - } - - - - template <> - inline - Vector<2>::symmetric_gradient_type - Vector<2>::symmetric_gradient (const unsigned int shape_function, - const unsigned int q_point) const - { - // this function works like in the case - // above - Assert (shape_function < fe_values.fe->dofs_per_cell, - ExcIndexRange (shape_function, 0, fe_values.fe->dofs_per_cell)); - Assert (fe_values.update_flags & update_gradients, - FEValuesBase<2>::ExcAccessToUninitializedField()); - - if (fe_values.fe->is_primitive() || - fe_values.fe->is_primitive(shape_function)) - { - const unsigned int - nonzero_component - = fe_values.fe->system_to_component_index(shape_function).first; - - // first get the one component of - // the nonsymmetrized gradient - // that is not zero - const Tensor<1,2> &grad - = fe_values.shape_gradients[fe_values. - shape_function_to_row_table[shape_function]][q_point]; - - // then form a symmetric tensor - // out of it. note that access to - // individual elements of a - // SymmetricTensor object is - // fairly slow. if we implemented - // the following in the naive - // way, we would therefore incur - // a very significant penalty: a - // preliminary version of the - // Stokes tutorial program would - // slow down from 17 to 23 - // seconds because of this single - // function! - // - // consequently, we try to be a - // bit smarter by already laying - // out the data in the right - // format and creating a - // symmetric tensor of it - switch (nonzero_component - first_vector_component) - { - case 0: - { - const double array[symmetric_gradient_type::n_independent_components] - = { grad[0], 0, grad[1]/2 }; - return symmetric_gradient_type(array); - } - - case 1: - { - const double array[symmetric_gradient_type::n_independent_components] - = { 0, grad[1], grad[0]/2 }; - return symmetric_gradient_type(array); - } - - default: - // not a shape - // function that - // shared in the - // components of - // the selected - // vector - return symmetric_gradient_type(); - } - } - else - { - return symmetrize (gradient(shape_function, q_point)); - } - } - - - - template <> + template inline - Vector<3>::symmetric_gradient_type - Vector<3>::symmetric_gradient (const unsigned int shape_function, - const unsigned int q_point) const + typename Vector::symmetric_gradient_type + Vector::symmetric_gradient (const unsigned int shape_function, + const unsigned int q_point) const { - // this function works like in the case - // above - Assert (shape_function < fe_values.fe->dofs_per_cell, - ExcIndexRange (shape_function, 0, fe_values.fe->dofs_per_cell)); - Assert (fe_values.update_flags & update_gradients, - FEValuesBase<3>::ExcAccessToUninitializedField()); - - if (fe_values.fe->is_primitive() || - fe_values.fe->is_primitive(shape_function)) - { - const unsigned int - nonzero_component - = fe_values.fe->system_to_component_index(shape_function).first; - - // first get the one component of - // the nonsymmetrized gradient - // that is not zero - const Tensor<1,3> &grad - = fe_values.shape_gradients[fe_values. - shape_function_to_row_table[shape_function]][q_point]; - - // then form a symmetric - // tensor out of it. note - // that access to individual - // elements of a - // SymmetricTensor object is - // fairly slow. if we - // implemented the following - // in the naive way, we would - // therefore incur a very - // significant penalty: a - // preliminary version of the - // Stokes tutorial program - // would slow down from 17 to - // 23 seconds because of this - // single function! - // - // consequently, we try to be - // a bit smarter by already - // laying out the data in the - // right format and creating - // a symmetric tensor of it - switch (nonzero_component-first_vector_component) - { - case 0: - { - const double array[symmetric_gradient_type::n_independent_components] - = { grad[0], 0, 0, grad[1]/2 , grad[2]/2, 0}; - return symmetric_gradient_type(array); - } - - case 1: - { - const double array[symmetric_gradient_type::n_independent_components] - = { 0, grad[1], 0, grad[0]/2, 0, grad[2]/2 }; - return symmetric_gradient_type(array); - } - - case 2: - { - const double array[symmetric_gradient_type::n_independent_components] - = { 0, 0, grad[2], 0, grad[0]/2, grad[1]/2 }; - return symmetric_gradient_type(array); - } - - default: - // not a shape - // function that - // shared in the - // components of - // the selected - // vector - return symmetric_gradient_type(); - } - } - else - { - return symmetrize (gradient(shape_function, q_point)); - } + return symmetrize (gradient(shape_function, q_point)); } } @@ -3922,23 +3622,31 @@ namespace FEValuesViews template inline -FEValuesViews::Scalar +const FEValuesViews::Scalar & FEValuesBase:: operator[] (const FEValuesExtractors::Scalar &scalar) const { - return FEValuesViews::Scalar (*this, scalar.component); + Assert (scalar.component < fe_values_views_cache.scalars.size(), + ExcIndexRange (scalar.component, + 0, fe_values_views_cache.scalars.size())); + + return fe_values_views_cache.scalars[scalar.component]; } template inline -FEValuesViews::Vector +const FEValuesViews::Vector & FEValuesBase:: operator[] (const FEValuesExtractors::Vector &vector) const { - return - FEValuesViews::Vector (*this, vector.first_vector_component); + Assert (vector.first_vector_component < + fe_values_views_cache.vectors.size(), + ExcIndexRange (vector.first_vector_component, + 0, fe_values_views_cache.vectors.size())); + + return fe_values_views_cache.vectors[vector.first_vector_component]; } diff --git a/deal.II/deal.II/source/fe/fe_values.cc b/deal.II/deal.II/source/fe/fe_values.cc index 808b78ca7f..6bee51751b 100644 --- a/deal.II/deal.II/source/fe/fe_values.cc +++ b/deal.II/deal.II/source/fe/fe_values.cc @@ -31,6 +31,217 @@ DEAL_II_NAMESPACE_OPEN + +namespace +{ + template + inline + std::vector + make_shape_function_to_row_table (const FiniteElement &fe) + { + std::vector shape_function_to_row_table (fe.dofs_per_cell); + unsigned int row = 0; + for (unsigned int i=0; i + Scalar::Scalar (const FEValuesBase &fe_values, + const unsigned int component) + : + fe_values (fe_values), + component (component), + is_nonzero_shape_function_component (fe_values.fe->dofs_per_cell), + row_index (fe_values.fe->dofs_per_cell) + { + Assert (component < fe_values.fe->n_components(), + ExcIndexRange(component, 0, fe_values.fe->n_components())); + + const std::vector shape_function_to_row_table + = make_shape_function_to_row_table (*fe_values.fe); + + for (unsigned int i=0; idofs_per_cell; ++i) + { + const bool is_primitive = (fe_values.fe->is_primitive() || + fe_values.fe->is_primitive(i)); + + if (is_primitive == true) + is_nonzero_shape_function_component[i] + = (component == + fe_values.fe->system_to_component_index(i).first); + else + is_nonzero_shape_function_component[i] + = (fe_values.fe->get_nonzero_components(i)[component] + == true); + + if (is_nonzero_shape_function_component[i] == true) + { + if (is_primitive == true) + row_index[i] = shape_function_to_row_table[i]; + else + row_index[i] = (shape_function_to_row_table[i] + + + std::count (fe_values.fe->get_nonzero_components(i).begin(), + fe_values.fe->get_nonzero_components(i).begin()+ + component, + true)); + } + else + row_index[i] = numbers::invalid_unsigned_int; + } + } + + + + template + Scalar::Scalar () + : + fe_values (*static_cast*>(0)), + component (numbers::invalid_unsigned_int) + {} + + + template + Scalar & + Scalar::operator= (const Scalar &) + { + // we shouldn't be copying these objects + Assert (false, ExcInternalError()); + return *this; + } + + + + template + Vector::Vector (const FEValuesBase &fe_values, + const unsigned int first_vector_component) + : + fe_values (fe_values), + first_vector_component (first_vector_component), + is_nonzero_shape_function_component (fe_values.fe->dofs_per_cell, + dim), + row_index (fe_values.fe->dofs_per_cell, + dim) + { + Assert (first_vector_component+dim-1 < fe_values.fe->n_components(), + ExcIndexRange(first_vector_component+dim-1, 0, + fe_values.fe->n_components())); + + const std::vector shape_function_to_row_table + = make_shape_function_to_row_table (*fe_values.fe); + + for (unsigned int d=0; ddofs_per_cell; ++i) + { + const bool is_primitive = (fe_values.fe->is_primitive() || + fe_values.fe->is_primitive(i)); + + if (is_primitive == true) + is_nonzero_shape_function_component[i][d] + = (component == + fe_values.fe->system_to_component_index(i).first); + else + is_nonzero_shape_function_component[i][d] + = (fe_values.fe->get_nonzero_components(i)[component] + == true); + + if (is_nonzero_shape_function_component[i][d] == true) + { + if (is_primitive == true) + row_index[i][d] = shape_function_to_row_table[i]; + else + row_index[i][d] = (shape_function_to_row_table[i] + + + std::count (fe_values.fe->get_nonzero_components(i).begin(), + fe_values.fe->get_nonzero_components(i).begin()+ + component, + true)); + } + else + row_index[i][d] = numbers::invalid_unsigned_int; + } + } + } + + + template + Vector::Vector () + : + fe_values (*static_cast*>(0)), + first_vector_component (numbers::invalid_unsigned_int) + {} + + + + template + Vector & + Vector::operator= (const Vector &) + { + // we shouldn't be copying these objects + Assert (false, ExcInternalError()); + return *this; + } +} + + +namespace internal +{ + namespace FEValuesViews + { + template + Cache::Cache (const FEValuesBase &fe_values) + { + const FiniteElement &fe = fe_values.get_fe(); + + // create the views objects. allocate a + // bunch of default-constructed ones + // then destroy them again and do + // in-place construction of those we + // actually want to use (copying stuff + // is wasteful and we can't do that + // anyway because the class has + // reference members) + const unsigned int n_scalars = fe.n_components(); + scalars.resize (n_scalars); + for (unsigned int component=0; component::~Scalar (); + new (&scalars[component]) + dealii::FEValuesViews::Scalar(fe_values, + component); + } + + const unsigned int n_vectors = (fe.n_components() >= dim ? + fe.n_components()-dim+1 : + 0); + vectors.resize (n_vectors); + for (unsigned int component=0; component::~Vector (); + new (&vectors[component]) + dealii::FEValuesViews::Vector(fe_values, + component); + } + } + } +} + + /* ---------------- FEValuesBase::CellIteratorBase --------- */ template @@ -407,19 +618,16 @@ FEValuesData::initialize (const unsigned int n_quadrature_p // from shape function number to // the rows in the tables denoting // its first non-zero - // component. with this also count - // the total number of non-zero - // components accumulated over all - // shape functions - this->shape_function_to_row_table.resize (fe.dofs_per_cell); - unsigned int row = 0; + // component + this->shape_function_to_row_table + = make_shape_function_to_row_table (fe); + + // count the total number of non-zero + // components accumulated over all shape + // functions + unsigned int n_nonzero_shape_components = 0; for (unsigned int i=0; ishape_function_to_row_table[i] = row; - row += fe.n_nonzero_components (i); - }; - - const unsigned int n_nonzero_shape_components = row; + n_nonzero_shape_components += fe.n_nonzero_components (i); Assert (n_nonzero_shape_components >= fe.dofs_per_cell, ExcInternalError()); @@ -483,7 +691,8 @@ FEValuesBase::FEValuesBase (const unsigned int n_q_points, mapping(&mapping), fe(&fe), mapping_data(0), - fe_data(0) + fe_data(0), + fe_values_views_cache (*this) { this->update_flags = flags; } @@ -1514,11 +1723,11 @@ FEValues::FEValues (const Mapping &mapping, const UpdateFlags update_flags) : FEValuesBase (q.size(), - fe.dofs_per_cell, - update_default, - mapping, - fe), - quadrature (q) + fe.dofs_per_cell, + update_default, + mapping, + fe), + quadrature (q) { initialize (update_flags); } -- 2.39.5