From d1da5b455ffa7ecb33848aca4c9e601a65351423 Mon Sep 17 00:00:00 2001 From: kronbichler Date: Fri, 14 Feb 2014 16:58:41 +0000 Subject: [PATCH] Remove template argument from FEEvaluationBase class. It is not needed there, only a derived class needs it. This is the first step to make FEEvaluation collaborate with some data generated by FEValues git-svn-id: https://svn.dealii.org/trunk@32489 0785d39b-7218-0410-832d-ea1e28bc413d --- .../deal.II/matrix_free/fe_evaluation.h | 1113 +++++++++-------- 1 file changed, 611 insertions(+), 502 deletions(-) diff --git a/deal.II/include/deal.II/matrix_free/fe_evaluation.h b/deal.II/include/deal.II/matrix_free/fe_evaluation.h index cc564b72b3..38ddc1bd49 100644 --- a/deal.II/include/deal.II/matrix_free/fe_evaluation.h +++ b/deal.II/include/deal.II/matrix_free/fe_evaluation.h @@ -39,14 +39,23 @@ namespace parallel +// forward declarations namespace internal { DeclException0 (ExcAccessToUninitializedField); template - void do_evaluate (FEEval &, const bool, const bool, const bool); + void do_evaluate (FEEval &, const bool, const bool, const bool, int2type<1>); template - void do_integrate (FEEval &, const bool, const bool); + void do_evaluate (FEEval &, const bool, const bool, const bool, int2type<2>); + template + void do_evaluate (FEEval &, const bool, const bool, const bool, int2type<3>); + template + void do_integrate (FEEval &, const bool, const bool, int2type<1>); + template + void do_integrate (FEEval &, const bool, const bool, int2type<2>); + template + void do_integrate (FEEval &, const bool, const bool, int2type<3>); } @@ -82,8 +91,7 @@ namespace internal * * @author Katharina Kormann and Martin Kronbichler, 2010, 2011 */ -template +template class FEEvaluationBase { public: @@ -92,8 +100,6 @@ public: typedef Tensor<1,n_components_,Tensor<1,dim,VectorizedArray > > gradient_type; static const unsigned int dimension = dim; static const unsigned int n_components = n_components_; - static const unsigned int dofs_per_cell = dofs_per_cell_; - static const unsigned int n_q_points = n_q_points_; /** * @name 1: General operations @@ -516,8 +522,10 @@ protected: * select the appropriate components. */ FEEvaluationBase (const MatrixFree &matrix_free, - const unsigned int fe_no = 0, - const unsigned int quad_no = 0); + const unsigned int fe_no, + const unsigned int quad_no, + const unsigned int dofs_per_cell, + const unsigned int n_q_points); /** * A unified function to read from and write into vectors based on the given @@ -551,7 +559,7 @@ protected: * or before distributing them into a result vector). The methods * get_dof_value() and submit_dof_value() read from or write to this field. */ - VectorizedArray values_dofs[n_components][dofs_per_cell>0?dofs_per_cell:1]; + VectorizedArray* values_dofs[n_components]; /** * This field stores the values of the finite element function on quadrature @@ -559,7 +567,7 @@ protected: * integrating. The methods get_value() and submit_value() access this * field. */ - VectorizedArray values_quad[n_components][n_q_points>0?n_q_points:1]; + VectorizedArray* values_quad[n_components]; /** * This field stores the gradients of the finite element function on @@ -568,14 +576,14 @@ protected: * some specializations like get_symmetric_gradient() or get_divergence()) * access this field. */ - VectorizedArray gradients_quad[n_components][dim][n_q_points>0?n_q_points:1]; + VectorizedArray* gradients_quad[n_components][dim]; /** * This field stores the Hessians of the finite element function on * quadrature points after applying unit cell transformations. The methods * get_hessian(), get_laplacian(), get_hessian_diagonal() access this field. */ - VectorizedArray hessians_quad[n_components][(dim*(dim+1))/2][n_q_points>0?n_q_points:1]; + VectorizedArray* hessians_quad[n_components][(dim*(dim+1))/2]; /** * Stores the number of the quadrature formula of the present cell. @@ -744,10 +752,8 @@ protected: * * @author Katharina Kormann and Martin Kronbichler, 2010, 2011 */ -template -class FEEvaluationAccess : - public FEEvaluationBase +template +class FEEvaluationAccess : public FEEvaluationBase { public: typedef Number number_type; @@ -755,10 +761,7 @@ public: typedef Tensor<1,n_components_,Tensor<1,dim,VectorizedArray > > gradient_type; static const unsigned int dimension = dim; static const unsigned int n_components = n_components_; - static const unsigned int dofs_per_cell = dofs_per_cell_; - static const unsigned int n_q_points = n_q_points_; - typedef FEEvaluationBase BaseClass; + typedef FEEvaluationBase BaseClass; protected: /** @@ -769,8 +772,10 @@ protected: * the appropriate components. */ FEEvaluationAccess (const MatrixFree &matrix_free, - const unsigned int fe_no = 0, - const unsigned int quad_no = 0); + const unsigned int fe_no, + const unsigned int quad_no, + const unsigned int dofs_per_cell, + const unsigned int n_q_points); }; @@ -784,18 +789,15 @@ protected: * * @author Katharina Kormann and Martin Kronbichler, 2010, 2011 */ -template -class FEEvaluationAccess : - public FEEvaluationBase +template +class FEEvaluationAccess : public FEEvaluationBase { public: typedef Number number_type; typedef VectorizedArray value_type; typedef Tensor<1,dim,VectorizedArray > gradient_type; static const unsigned int dimension = dim; - static const unsigned int dofs_per_cell = dofs_per_cell_; - static const unsigned int n_q_points = n_q_points_; - typedef FEEvaluationBase BaseClass; + typedef FEEvaluationBase BaseClass; /** * Returns the value stored for the local degree of freedom with index @p @@ -892,8 +894,10 @@ protected: * the appropriate components. */ FEEvaluationAccess (const MatrixFree &matrix_free, - const unsigned int fe_no = 0, - const unsigned int quad_no = 0); + const unsigned int fe_no, + const unsigned int quad_no, + const unsigned int dofs_per_cell, + const unsigned int n_q_points); }; @@ -907,9 +911,8 @@ protected: * * @author Katharina Kormann and Martin Kronbichler, 2010, 2011 */ -template -class FEEvaluationAccess : - public FEEvaluationBase +template +class FEEvaluationAccess : public FEEvaluationBase { public: typedef Number number_type; @@ -917,9 +920,7 @@ public: typedef Tensor<2,dim,VectorizedArray > gradient_type; static const unsigned int dimension = dim; static const unsigned int n_components = dim; - static const unsigned int dofs_per_cell = dofs_per_cell_; - static const unsigned int n_q_points = n_q_points_; - typedef FEEvaluationBase BaseClass; + typedef FEEvaluationBase BaseClass; /** * Returns the gradient of a finite element function at quadrature point @@ -1024,8 +1025,10 @@ protected: * the appropriate components. */ FEEvaluationAccess (const MatrixFree &matrix_free, - const unsigned int fe_no = 0, - const unsigned int quad_no = 0); + const unsigned int fe_no, + const unsigned int quad_no, + const unsigned int dofs_per_cell, + const unsigned int n_q_points); }; @@ -1071,24 +1074,17 @@ protected: */ template -class FEEvaluationGeneral : - public FEEvaluationAccess::value, - Utilities::fixed_int_power::value, - n_components_,Number> +class FEEvaluationGeneral : public FEEvaluationAccess { public: - typedef FEEvaluationAccess::value, - Utilities::fixed_int_power::value, - n_components_, Number> BaseClass; + typedef FEEvaluationAccess BaseClass; typedef Number number_type; typedef typename BaseClass::value_type value_type; typedef typename BaseClass::gradient_type gradient_type; static const unsigned int dimension = dim; static const unsigned int n_components = n_components_; - static const unsigned int dofs_per_cell = BaseClass::dofs_per_cell; - static const unsigned int n_q_points = BaseClass::n_q_points; + static const unsigned int dofs_per_cell = Utilities::fixed_int_power::value; + static const unsigned int n_q_points = Utilities::fixed_int_power::value; /** * Constructor. Takes all data stored in MatrixFree. If applied to problems @@ -1170,9 +1166,14 @@ protected: * Friend declaration. */ template friend void - internal::do_evaluate (FEEval &, const bool, const bool, const bool); + internal::do_evaluate (FEEval &, const bool, const bool, const bool, internal::int2type); template friend void - internal::do_integrate (FEEval &, const bool, const bool); + internal::do_integrate (FEEval &, const bool, const bool, internal::int2type); + + /** + * Internally stored variables for the different data fields. + */ + VectorizedArray my_data_array[n_components*(dofs_per_cell+(dim*dim+2*dim+1)*n_q_points)]; }; @@ -1312,9 +1313,9 @@ protected: * Friend declarations. */ template friend void - internal::do_evaluate (FEEval &, const bool, const bool, const bool); + internal::do_evaluate (FEEval &, const bool, const bool, const bool, internal::int2type); template friend void - internal::do_integrate (FEEval &, const bool, const bool); + internal::do_integrate (FEEval &, const bool, const bool, internal::int2type); }; @@ -1426,21 +1427,22 @@ protected: /*----------------------- FEEvaluationBase ----------------------------------*/ -template +template inline -FEEvaluationBase +FEEvaluationBase ::FEEvaluationBase (const MatrixFree &data_in, const unsigned int fe_no_in, - const unsigned int quad_no_in) + const unsigned int quad_no_in, + const unsigned int dofs_per_cell, + const unsigned int n_q_points) : quad_no (quad_no_in), n_fe_components (data_in.get_dof_info(fe_no_in).n_components), active_fe_index (data_in.get_dof_info(fe_no_in).fe_index_from_dofs_per_cell - (dofs_per_cell_ * n_fe_components)), + (dofs_per_cell * n_fe_components)), active_quad_index (data_in.get_mapping_info(). mapping_data_gen[quad_no_in]. - quad_index_from_n_q_points(n_q_points_)), + quad_index_from_n_q_points(n_q_points)), matrix_info (data_in), dof_info (data_in.get_dof_info(fe_no_in)), mapping_info (data_in.get_mapping_info()), @@ -1459,10 +1461,23 @@ FEEvaluationBase cell_type (internal::MatrixFreeFunctions::undefined), cell_data_number (0) { + for (unsigned int c=0; c::n_array_elements); + AssertDimension (data.dofs_per_cell, + dof_info.dofs_per_cell[active_fe_index]/n_fe_components); + AssertDimension (data.n_q_points, + mapping_info.mapping_data_gen[quad_no].n_q_points[active_quad_index]); Assert (n_fe_components == 1 || n_components == 1 || n_components == n_fe_components, @@ -1477,11 +1492,10 @@ FEEvaluationBase -template +template inline void -FEEvaluationBase +FEEvaluationBase ::reinit (const unsigned int cell_in) { AssertIndexRange (cell_in, dof_info.row_starts.size()-1); @@ -1552,11 +1566,10 @@ FEEvaluationBase -template +template inline unsigned int -FEEvaluationBase +FEEvaluationBase ::get_cell_data_number () const { Assert (cell != numbers::invalid_unsigned_int, ExcNotInitialized()); @@ -1565,11 +1578,10 @@ FEEvaluationBase -template +template inline internal::MatrixFreeFunctions::CellType -FEEvaluationBase +FEEvaluationBase ::get_cell_type () const { Assert (cell != numbers::invalid_unsigned_int, ExcNotInitialized()); @@ -1801,13 +1813,11 @@ namespace internal - -template +template template inline void -FEEvaluationBase +FEEvaluationBase ::read_write_operation (const VectorOperation &operation, VectorType *src[]) const { @@ -1831,6 +1841,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 n_irreg_components_filled = dof_info.row_starts[cell][2]; const bool at_irregular_cell = n_irreg_components_filled > 0; @@ -2123,12 +2134,11 @@ FEEvaluationBase -template +template template inline void -FEEvaluationBase +FEEvaluationBase ::read_dof_values (const VectorType &src) { // select between block vectors and non-block vectors. Note that the number @@ -2148,12 +2158,11 @@ FEEvaluationBase -template +template template inline void -FEEvaluationBase +FEEvaluationBase ::read_dof_values (const std::vector &src, const unsigned int first_index) { @@ -2177,12 +2186,11 @@ FEEvaluationBase -template +template template inline void -FEEvaluationBase +FEEvaluationBase ::read_dof_values (const std::vector &src, const unsigned int first_index) { @@ -2206,12 +2214,11 @@ FEEvaluationBase -template +template template inline void -FEEvaluationBase +FEEvaluationBase ::read_dof_values_plain (const VectorType &src) { // select between block vectors and non-block vectors. Note that the number @@ -2226,12 +2233,11 @@ FEEvaluationBase -template +template template inline void -FEEvaluationBase +FEEvaluationBase ::read_dof_values_plain (const std::vector &src, const unsigned int first_index) { @@ -2248,12 +2254,11 @@ FEEvaluationBase -template +template template inline void -FEEvaluationBase +FEEvaluationBase ::read_dof_values_plain (const std::vector &src, const unsigned int first_index) { @@ -2270,12 +2275,11 @@ FEEvaluationBase -template +template template inline void -FEEvaluationBase +FEEvaluationBase ::distribute_local_to_global (VectorType &dst) const { Assert (dof_values_initialized==true, @@ -2294,12 +2298,11 @@ FEEvaluationBase -template +template template inline void -FEEvaluationBase +FEEvaluationBase ::distribute_local_to_global (std::vector &dst, const unsigned int first_index) const { @@ -2321,12 +2324,11 @@ FEEvaluationBase -template +template template inline void -FEEvaluationBase +FEEvaluationBase ::distribute_local_to_global (std::vector &dst, const unsigned int first_index) const { @@ -2348,12 +2350,11 @@ FEEvaluationBase -template +template template inline void -FEEvaluationBase +FEEvaluationBase ::set_dof_values (VectorType &dst) const { Assert (dof_values_initialized==true, @@ -2372,12 +2373,11 @@ FEEvaluationBase -template +template template inline void -FEEvaluationBase +FEEvaluationBase ::set_dof_values (std::vector &dst, const unsigned int first_index) const { @@ -2400,12 +2400,11 @@ FEEvaluationBase -template +template template inline void -FEEvaluationBase +FEEvaluationBase ::set_dof_values (std::vector &dst, const unsigned int first_index) const { @@ -2428,12 +2427,11 @@ FEEvaluationBase -template +template template inline void -FEEvaluationBase +FEEvaluationBase ::read_dof_values_plain (const VectorType *src[]) { // this is different from the other three operations because we do not use @@ -2447,6 +2445,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 n_irreg_components_filled = dof_info.row_starts[cell][2]; const bool at_irregular_cell = n_irreg_components_filled > 0; @@ -2548,11 +2547,10 @@ FEEvaluationBase /*------------------------------ access to data fields ----------------------*/ -template +template inline const VectorizedArray * -FEEvaluationBase:: +FEEvaluationBase:: begin_dof_values () const { return &values_dofs[0][0]; @@ -2560,11 +2558,10 @@ begin_dof_values () const -template +template inline VectorizedArray * -FEEvaluationBase:: +FEEvaluationBase:: begin_dof_values () { #ifdef DEBUG @@ -2575,11 +2572,10 @@ begin_dof_values () -template +template inline const VectorizedArray * -FEEvaluationBase:: +FEEvaluationBase:: begin_values () const { Assert (values_quad_initialized || values_quad_submitted, @@ -2589,11 +2585,10 @@ begin_values () const -template +template inline VectorizedArray * -FEEvaluationBase:: +FEEvaluationBase:: begin_values () { #ifdef DEBUG @@ -2604,11 +2599,10 @@ begin_values () -template +template inline const VectorizedArray * -FEEvaluationBase:: +FEEvaluationBase:: begin_gradients () const { Assert (gradients_quad_initialized || gradients_quad_submitted, @@ -2618,11 +2612,10 @@ begin_gradients () const -template +template inline VectorizedArray * -FEEvaluationBase:: +FEEvaluationBase:: begin_gradients () { #ifdef DEBUG @@ -2633,11 +2626,10 @@ begin_gradients () -template +template inline const VectorizedArray * -FEEvaluationBase:: +FEEvaluationBase:: begin_hessians () const { Assert (hessians_quad_initialized, ExcNotInitialized()); @@ -2646,11 +2638,10 @@ begin_hessians () const -template +template inline VectorizedArray * -FEEvaluationBase:: +FEEvaluationBase:: begin_hessians () { return &hessians_quad[0][0][0]; @@ -2658,14 +2649,13 @@ begin_hessians () -template +template inline Tensor<1,n_components_,VectorizedArray > -FEEvaluationBase +FEEvaluationBase ::get_dof_value (const unsigned int dof) const { - AssertIndexRange (dof, dofs_per_cell); + AssertIndexRange (dof, this->data.dofs_per_cell); Tensor<1,n_components_,VectorizedArray > return_value (false); for (unsigned int comp=0; compvalues_dofs[comp][dof]; @@ -2674,16 +2664,15 @@ FEEvaluationBase -template +template inline Tensor<1,n_components_,VectorizedArray > -FEEvaluationBase +FEEvaluationBase ::get_value (const unsigned int q_point) const { Assert (this->values_quad_initialized==true, internal::ExcAccessToUninitializedField()); - AssertIndexRange (q_point, n_q_points); + AssertIndexRange (q_point, this->data.n_q_points); Tensor<1,n_components_,VectorizedArray > return_value (false); for (unsigned int comp=0; compvalues_quad[comp][q_point]; @@ -2692,16 +2681,15 @@ FEEvaluationBase -template +template inline Tensor<1,n_components_,Tensor<1,dim,VectorizedArray > > -FEEvaluationBase +FEEvaluationBase ::get_gradient (const unsigned int q_point) const { Assert (this->gradients_quad_initialized==true, internal::ExcAccessToUninitializedField()); - AssertIndexRange (q_point, n_q_points); + AssertIndexRange (q_point, this->data.n_q_points); Tensor<1,n_components_,Tensor<1,dim,VectorizedArray > > grad_out (false); @@ -2740,57 +2728,68 @@ namespace internal { // compute tmp = hess_unit(u) * J^T. do this manually because we do not // store the lower diagonal because of symmetry - template + template inline void - hessian_unit_times_jac (const Tensor<2,dim,VectorizedArray > &jac, - const VectorizedArray hessians_quad[][n_q_points], - const unsigned int q_point, - VectorizedArray tmp[dim][dim]) + hessian_unit_times_jac (const Tensor<2,1,VectorizedArray > &jac, + const VectorizedArray *const hessians_quad[1], + const unsigned int q_point, + VectorizedArray (&tmp)[1][1]) { - for (unsigned int d=0; d + inline + void + hessian_unit_times_jac (const Tensor<2,2,VectorizedArray > &jac, + const VectorizedArray *const hessians_quad[3], + const unsigned int q_point, + VectorizedArray (&tmp)[2][2]) + { + for (unsigned int d=0; d<2; ++d) { - switch (dim) - { - case 1: - tmp[0][0] = jac[0][0] * hessians_quad[0][q_point]; - break; - case 2: - tmp[0][d] = (jac[d][0] * hessians_quad[0][q_point] + - jac[d][1] * hessians_quad[2][q_point]); - tmp[1][d] = (jac[d][0] * hessians_quad[2][q_point] + - jac[d][1] * hessians_quad[1][q_point]); - break; - case 3: - tmp[0][d] = (jac[d][0] * hessians_quad[0][q_point] + - jac[d][1] * hessians_quad[3][q_point] + - jac[d][2] * hessians_quad[4][q_point]); - tmp[1][d] = (jac[d][0] * hessians_quad[3][q_point] + - jac[d][1] * hessians_quad[1][q_point] + - jac[d][2] * hessians_quad[5][q_point]); - tmp[2][d] = (jac[d][0] * hessians_quad[4][q_point] + - jac[d][1] * hessians_quad[5][q_point] + - jac[d][2] * hessians_quad[2][q_point]); - break; - default: - Assert (false, ExcNotImplemented()); - } + tmp[0][d] = (jac[d][0] * hessians_quad[0][q_point] + + jac[d][1] * hessians_quad[2][q_point]); + tmp[1][d] = (jac[d][0] * hessians_quad[2][q_point] + + jac[d][1] * hessians_quad[1][q_point]); + } + } + + template + inline + void + hessian_unit_times_jac (const Tensor<2,3,VectorizedArray > &jac, + const VectorizedArray *const hessians_quad[6], + const unsigned int q_point, + VectorizedArray (&tmp)[3][3]) + { + for (unsigned int d=0; d<3; ++d) + { + tmp[0][d] = (jac[d][0] * hessians_quad[0][q_point] + + jac[d][1] * hessians_quad[3][q_point] + + jac[d][2] * hessians_quad[4][q_point]); + tmp[1][d] = (jac[d][0] * hessians_quad[3][q_point] + + jac[d][1] * hessians_quad[1][q_point] + + jac[d][2] * hessians_quad[5][q_point]); + tmp[2][d] = (jac[d][0] * hessians_quad[4][q_point] + + jac[d][1] * hessians_quad[5][q_point] + + jac[d][2] * hessians_quad[2][q_point]); } } } -template +template inline Tensor<1,n_components_,Tensor<2,dim,VectorizedArray > > -FEEvaluationBase +FEEvaluationBase ::get_hessian (const unsigned int q_point) const { Assert (this->hessians_quad_initialized==true, internal::ExcAccessToUninitializedField()); - AssertIndexRange (q_point, n_q_points); + AssertIndexRange (q_point, this->data.n_q_points); Tensor<2,dim,VectorizedArray > hessian_out [n_components]; @@ -2907,16 +2906,15 @@ FEEvaluationBase -template +template inline Tensor<1,n_components_,Tensor<1,dim,VectorizedArray > > -FEEvaluationBase +FEEvaluationBase ::get_hessian_diagonal (const unsigned int q_point) const { Assert (this->hessians_quad_initialized==true, internal::ExcAccessToUninitializedField()); - AssertIndexRange (q_point, n_q_points); + AssertIndexRange (q_point, this->data.n_q_points); Tensor<1,n_components_,Tensor<1,dim,VectorizedArray > > hessian_out (false); @@ -2986,16 +2984,15 @@ FEEvaluationBase -template +template inline Tensor<1,n_components_,VectorizedArray > -FEEvaluationBase +FEEvaluationBase ::get_laplacian (const unsigned int q_point) const { Assert (this->hessians_quad_initialized==true, internal::ExcAccessToUninitializedField()); - AssertIndexRange (q_point, n_q_points); + AssertIndexRange (q_point, this->data.n_q_points); Tensor<1,n_components_,VectorizedArray > laplacian_out (false); const Tensor<1,n_components_,Tensor<1,dim,VectorizedArray > > hess_diag = get_hessian_diagonal(q_point); @@ -3010,35 +3007,33 @@ FEEvaluationBase -template +template inline void -FEEvaluationBase +FEEvaluationBase ::submit_dof_value (const Tensor<1,n_components_,VectorizedArray > val_in, const unsigned int dof) { #ifdef DEBUG this->dof_values_initialized = true; #endif - AssertIndexRange (dof, dofs_per_cell); + AssertIndexRange (dof, this->data.dofs_per_cell); for (unsigned int comp=0; compvalues_dofs[comp][dof] = val_in[comp]; } -template +template inline void -FEEvaluationBase +FEEvaluationBase ::submit_value (const Tensor<1,n_components_,VectorizedArray > val_in, const unsigned int q_point) { #ifdef DEBUG Assert (this->cell != numbers::invalid_unsigned_int, ExcNotInitialized()); - AssertIndexRange (q_point, n_q_points); + AssertIndexRange (q_point, this->data.n_q_points); this->values_quad_submitted = true; #endif if (this->cell_type == internal::MatrixFreeFunctions::general) @@ -3057,18 +3052,17 @@ FEEvaluationBase -template +template inline void -FEEvaluationBase +FEEvaluationBase ::submit_gradient (const Tensor<1,n_components_, Tensor<1,dim,VectorizedArray > >grad_in, const unsigned int q_point) { #ifdef DEBUG Assert (this->cell != numbers::invalid_unsigned_int, ExcNotInitialized()); - AssertIndexRange (q_point, n_q_points); + AssertIndexRange (q_point, this->data.n_q_points); this->gradients_quad_submitted = true; #endif if (this->cell_type == internal::MatrixFreeFunctions::cartesian) @@ -3100,11 +3094,10 @@ FEEvaluationBase -template +template inline Tensor<1,n_components_,VectorizedArray > -FEEvaluationBase +FEEvaluationBase ::integrate_value () const { #ifdef DEBUG @@ -3115,6 +3108,7 @@ FEEvaluationBase Tensor<1,n_components_,VectorizedArray > return_value (false); for (unsigned int comp=0; compvalues_quad[comp][0]; + const unsigned int n_q_points = this->data.n_q_points; for (unsigned int q=1; qvalues_quad[comp][q]; @@ -3126,16 +3120,17 @@ FEEvaluationBase /*----------------------- FEEvaluationAccess --------------------------------*/ -template +template inline -FEEvaluationAccess +FEEvaluationAccess ::FEEvaluationAccess (const MatrixFree &data_in, const unsigned int fe_no, - const unsigned int quad_no_in) + const unsigned int quad_no_in, + const unsigned int dofs_per_cell, + const unsigned int n_q_points) : - FEEvaluationBase - (data_in, fe_no, quad_no_in) + FEEvaluationBase + (data_in, fe_no, quad_no_in, dofs_per_cell, n_q_points) {} @@ -3144,49 +3139,51 @@ FEEvaluationAccess /*-------------------- FEEvaluationAccess scalar ----------------------------*/ -template +template inline -FEEvaluationAccess +FEEvaluationAccess ::FEEvaluationAccess (const MatrixFree &data_in, const unsigned int fe_no, - const unsigned int quad_no_in) + const unsigned int quad_no_in, + const unsigned int dofs_per_cell, + const unsigned int n_q_points) : - FEEvaluationBase - (data_in, fe_no, quad_no_in) + FEEvaluationBase + (data_in, fe_no, quad_no_in, dofs_per_cell, n_q_points) {} -template +template inline VectorizedArray -FEEvaluationAccess +FEEvaluationAccess ::get_dof_value (const unsigned int dof) const { - AssertIndexRange (dof, dofs_per_cell); + AssertIndexRange (dof, this->data.dofs_per_cell); return this->values_dofs[0][dof]; } -template +template inline VectorizedArray -FEEvaluationAccess +FEEvaluationAccess ::get_value (const unsigned int q_point) const { Assert (this->values_quad_initialized==true, internal::ExcAccessToUninitializedField()); - AssertIndexRange (q_point, n_q_points); + AssertIndexRange (q_point, this->data.n_q_points); return this->values_quad[0][q_point]; } -template +template inline Tensor<1,dim,VectorizedArray > -FEEvaluationAccess +FEEvaluationAccess ::get_gradient (const unsigned int q_point) const { // could use the base class gradient, but that involves too many inefficient @@ -3194,7 +3191,7 @@ FEEvaluationAccess Assert (this->gradients_quad_initialized==true, internal::ExcAccessToUninitializedField()); - AssertIndexRange (q_point, n_q_points); + AssertIndexRange (q_point, this->data.n_q_points); Tensor<1,dim,VectorizedArray > grad_out (false); @@ -3223,10 +3220,10 @@ FEEvaluationAccess -template +template inline Tensor<2,dim,VectorizedArray > -FEEvaluationAccess +FEEvaluationAccess ::get_hessian (const unsigned int q_point) const { return BaseClass::get_hessian(q_point)[0]; @@ -3234,10 +3231,10 @@ FEEvaluationAccess -template +template inline Tensor<1,dim,VectorizedArray > -FEEvaluationAccess +FEEvaluationAccess ::get_hessian_diagonal (const unsigned int q_point) const { return BaseClass::get_hessian_diagonal(q_point)[0]; @@ -3245,10 +3242,10 @@ FEEvaluationAccess -template +template inline VectorizedArray -FEEvaluationAccess +FEEvaluationAccess ::get_laplacian (const unsigned int q_point) const { return BaseClass::get_laplacian(q_point)[0]; @@ -3256,32 +3253,32 @@ FEEvaluationAccess -template +template inline void -FEEvaluationAccess +FEEvaluationAccess ::submit_dof_value (const VectorizedArray val_in, const unsigned int dof) { #ifdef DEBUG this->dof_values_initialized = true; - AssertIndexRange (dof, dofs_per_cell); + AssertIndexRange (dof, this->data.dofs_per_cell); #endif this->values_dofs[0][dof] = val_in; } -template +template inline void -FEEvaluationAccess +FEEvaluationAccess ::submit_value (const VectorizedArray val_in, const unsigned int q_point) { #ifdef DEBUG Assert (this->cell != numbers::invalid_unsigned_int, ExcNotInitialized()); - AssertIndexRange (q_point, n_q_points); + AssertIndexRange (q_point, this->data.n_q_points); this->values_quad_submitted = true; #endif if (this->cell_type == internal::MatrixFreeFunctions::general) @@ -3298,16 +3295,16 @@ FEEvaluationAccess -template +template inline void -FEEvaluationAccess +FEEvaluationAccess ::submit_gradient (const Tensor<1,dim,VectorizedArray > grad_in, const unsigned int q_point) { #ifdef DEBUG Assert (this->cell != numbers::invalid_unsigned_int, ExcNotInitialized()); - AssertIndexRange (q_point, n_q_points); + AssertIndexRange (q_point, this->data.n_q_points); this->gradients_quad_submitted = true; #endif if (this->cell_type == internal::MatrixFreeFunctions::cartesian) @@ -3339,10 +3336,10 @@ FEEvaluationAccess -template +template inline VectorizedArray -FEEvaluationAccess +FEEvaluationAccess ::integrate_value () const { return BaseClass::integrate_value()[0]; @@ -3354,23 +3351,25 @@ FEEvaluationAccess /*----------------- FEEvaluationAccess vector-valued ------------------------*/ -template +template inline -FEEvaluationAccess +FEEvaluationAccess ::FEEvaluationAccess (const MatrixFree &data_in, const unsigned int fe_no, - const unsigned int quad_no_in) + const unsigned int quad_no_in, + const unsigned int dofs_per_cell, + const unsigned int n_q_points) : - FEEvaluationBase - (data_in, fe_no, quad_no_in) + FEEvaluationBase + (data_in, fe_no, quad_no_in, dofs_per_cell, n_q_points) {} -template +template inline Tensor<2,dim,VectorizedArray > -FEEvaluationAccess +FEEvaluationAccess ::get_gradient (const unsigned int q_point) const { return BaseClass::get_gradient (q_point); @@ -3378,15 +3377,15 @@ FEEvaluationAccess -template +template inline VectorizedArray -FEEvaluationAccess +FEEvaluationAccess ::get_divergence (const unsigned int q_point) const { Assert (this->gradients_quad_initialized==true, internal::ExcAccessToUninitializedField()); - AssertIndexRange (q_point, n_q_points); + AssertIndexRange (q_point, this->data.n_q_points); VectorizedArray divergence; @@ -3417,10 +3416,10 @@ FEEvaluationAccess -template +template inline SymmetricTensor<2,dim,VectorizedArray > -FEEvaluationAccess +FEEvaluationAccess ::get_symmetric_gradient (const unsigned int q_point) const { // copy from generic function into dim-specialization function @@ -3453,10 +3452,10 @@ FEEvaluationAccess -template +template inline Tensor<1,dim==2?1:dim,VectorizedArray > -FEEvaluationAccess +FEEvaluationAccess ::get_curl (const unsigned int q_point) const { // copy from generic function into dim-specialization function @@ -3484,39 +3483,39 @@ FEEvaluationAccess -template +template inline Tensor<2,dim,VectorizedArray > -FEEvaluationAccess +FEEvaluationAccess ::get_hessian_diagonal (const unsigned int q_point) const { Assert (this->hessians_quad_initialized==true, internal::ExcAccessToUninitializedField()); - AssertIndexRange (q_point, n_q_points); + AssertIndexRange (q_point, this->data.n_q_points); return BaseClass::get_hessian_diagonal (q_point); } -template +template inline Tensor<3,dim,VectorizedArray > -FEEvaluationAccess +FEEvaluationAccess ::get_hessian (const unsigned int q_point) const { Assert (this->hessians_quad_initialized==true, internal::ExcAccessToUninitializedField()); - AssertIndexRange (q_point, n_q_points); + AssertIndexRange (q_point, this->data.n_q_points); return BaseClass::get_hessian(q_point); } -template +template inline void -FEEvaluationAccess +FEEvaluationAccess ::submit_gradient (const Tensor<2,dim,VectorizedArray > grad_in, const unsigned int q_point) { @@ -3525,10 +3524,10 @@ FEEvaluationAccess -template +template inline void -FEEvaluationAccess +FEEvaluationAccess ::submit_gradient (const Tensor<1,dim,Tensor<1,dim,VectorizedArray > > grad_in, const unsigned int q_point) @@ -3536,17 +3535,18 @@ FEEvaluationAccess BaseClass::submit_gradient(grad_in, q_point); } -template + + +template inline void -FEEvaluationAccess +FEEvaluationAccess ::submit_divergence (const VectorizedArray div_in, const unsigned int q_point) { #ifdef DEBUG Assert (this->cell != numbers::invalid_unsigned_int, ExcNotInitialized()); - AssertIndexRange (q_point, n_q_points); + AssertIndexRange (q_point, this->data.n_q_points); this->gradients_quad_submitted = true; #endif if (this->cell_type == internal::MatrixFreeFunctions::cartesian) @@ -3583,10 +3583,10 @@ FEEvaluationAccess -template +template inline void -FEEvaluationAccess +FEEvaluationAccess ::submit_symmetric_gradient(const SymmetricTensor<2,dim,VectorizedArray > sym_grad, const unsigned int q_point) @@ -3596,7 +3596,7 @@ FEEvaluationAccess // that saves some operations #ifdef DEBUG Assert (this->cell != numbers::invalid_unsigned_int, ExcNotInitialized()); - AssertIndexRange (q_point, n_q_points); + AssertIndexRange (q_point, this->data.n_q_points); this->gradients_quad_submitted = true; #endif if (this->cell_type == internal::MatrixFreeFunctions::cartesian) @@ -3648,10 +3648,10 @@ FEEvaluationAccess -template +template inline void -FEEvaluationAccess +FEEvaluationAccess ::submit_curl (const Tensor<1,dim==2?1:dim,VectorizedArray > curl, const unsigned int q_point) { @@ -3692,8 +3692,22 @@ FEEvaluationGeneral const unsigned int fe_no, const unsigned int quad_no) : - BaseClass (data_in, fe_no, quad_no) + BaseClass (data_in, fe_no, quad_no, dofs_per_cell, n_q_points) { + // set the pointers to the correct position in the data array + for (unsigned int c=0; cvalues_dofs[c] = &my_data_array[c*dofs_per_cell]; + this->values_quad[c] = &my_data_array[n_components*dofs_per_cell+c*n_q_points]; + for (unsigned int d=0; dgradients_quad[c][d] = &my_data_array[n_components*(dofs_per_cell+n_q_points) + + + (c*dim+d)*n_q_points]; + for (unsigned int d=0; d<(dim*dim+dim)/2; ++d) + this->hessians_quad[c][d] = &my_data_array[n_components*((dim+1)*n_q_points+dofs_per_cell) + + + (c*(dim*dim+dim)+d)*n_q_points]; + } #ifdef DEBUG // print error message when the dimensions do not match. Propose a possible // fix @@ -4866,15 +4880,59 @@ namespace internal // This performs the evaluation of function values, gradients and Hessians // for tensor-product finite elements. The operation is used for both // FEEvaluationGeneral and FEEvaluation, which provide different functions - // apply_values, apply_gradients in the individual coordinate directions + // apply_values, apply_gradients in the individual coordinate directions. + // use different versions for 1d, 2d, 3d to avoid nasty compiler warnings template inline void do_evaluate (FEEval &fe_eval, const bool evaluate_val, const bool evaluate_grad, - const bool evaluate_lapl) + const bool evaluate_lapl, + internal::int2type<1>) { + AssertDimension(FEEval::dimension, 1); + Assert (fe_eval.cell != numbers::invalid_unsigned_int, + ExcNotInitialized()); + Assert (fe_eval.dof_values_initialized == true, + internal::ExcAccessToUninitializedField()); + + const unsigned int n_components = FEEval::n_components; + + for (unsigned int c=0; c + (fe_eval.values_dofs[c], fe_eval.values_quad[c]); + if (evaluate_grad == true) + fe_eval.template apply_gradients<0,true,false> + (fe_eval.values_dofs[c], fe_eval.gradients_quad[c][0]); + if (evaluate_lapl == true) + fe_eval.template apply_hessians<0,true,false> + (fe_eval.values_dofs[c], fe_eval.hessians_quad[c][0]); + } + +#ifdef DEBUG + if (evaluate_val == true) + fe_eval.values_quad_initialized = true; + if (evaluate_grad == true) + fe_eval.gradients_quad_initialized = true; + if (evaluate_lapl == true) + fe_eval.hessians_quad_initialized = true; +#endif + } + + + template + inline + void + do_evaluate (FEEval &fe_eval, + const bool evaluate_val, + const bool evaluate_grad, + const bool evaluate_lapl, + internal::int2type<2>) + { + AssertDimension(FEEval::dimension, 2); Assert (fe_eval.cell != numbers::invalid_unsigned_int, ExcNotInitialized()); Assert (fe_eval.dof_values_initialized == true, @@ -4883,163 +4941,171 @@ namespace internal const unsigned int temp_size = FEEval::dofs_per_cell > FEEval::n_q_points ? FEEval::dofs_per_cell : FEEval::n_q_points; const unsigned int n_components = FEEval::n_components; - const unsigned int dim = FEEval::dimension; for (unsigned int c=0; c temp1[temp_size]; VectorizedArray temp2[temp_size]; - switch (dim) + // grad x + if (evaluate_grad == true) { - case 3: - - if (evaluate_grad == true) - { - // grad x - fe_eval.template apply_gradients<0,true,false> + fe_eval.template apply_gradients<0,true,false> + (fe_eval.values_dofs[c], temp1); + fe_eval.template apply_values<1,true,false> + (temp1, fe_eval.gradients_quad[c][0]); + } + if (evaluate_lapl == true) + { + // grad xy + if (evaluate_grad == false) + fe_eval.template apply_gradients<0,true,false> (fe_eval.values_dofs[c], temp1); - fe_eval.template apply_values<1,true,false> - (temp1, temp2); - fe_eval.template apply_values<2,true,false> - (temp2, fe_eval.gradients_quad[c][0]); - } + fe_eval.template apply_gradients<1,true,false> + (temp1, fe_eval.hessians_quad[c][2]); - if (evaluate_lapl == true) - { - // grad xz - if (evaluate_grad == false) - { - fe_eval.template apply_gradients<0,true,false> - (fe_eval.values_dofs[c], temp1); - fe_eval.template apply_values<1,true,false> - (temp1, temp2); - } - fe_eval.template apply_gradients<2,true,false> - (temp2, fe_eval.hessians_quad[c][4]); + // grad xx + fe_eval.template apply_hessians<0,true,false> + (fe_eval.values_dofs[c], temp1); + fe_eval.template apply_values<1,true,false> + (temp1, fe_eval.hessians_quad[c][0]); + } - // grad xy - fe_eval.template apply_gradients<1,true,false> - (temp1, temp2); - fe_eval.template apply_values<2,true,false> - (temp2, fe_eval.hessians_quad[c][3]); + // grad y + fe_eval.template apply_values<0,true,false> + (fe_eval.values_dofs[c], temp1); + if (evaluate_grad == true) + fe_eval.template apply_gradients<1,true,false> + (temp1, fe_eval.gradients_quad[c][1]); + + // grad yy + if (evaluate_lapl == true) + fe_eval.template apply_hessians<1,true,false> + (temp1, fe_eval.hessians_quad[c][1]); + + // val: can use values applied in x + if (evaluate_val == true) + fe_eval.template apply_values<1,true,false> + (temp1, fe_eval.values_quad[c]); + } - // grad xx - fe_eval.template apply_hessians<0,true,false> - (fe_eval.values_dofs[c], temp1); - fe_eval.template apply_values<1,true,false> - (temp1, temp2); - fe_eval.template apply_values<2,true,false> - (temp2, fe_eval.hessians_quad[c][0]); - } +#ifdef DEBUG + if (evaluate_val == true) + fe_eval.values_quad_initialized = true; + if (evaluate_grad == true) + fe_eval.gradients_quad_initialized = true; + if (evaluate_lapl == true) + fe_eval.hessians_quad_initialized = true; +#endif + } - // grad y - fe_eval.template apply_values<0,true,false> - (fe_eval.values_dofs[c], temp1); - if (evaluate_grad == true) - { - fe_eval.template apply_gradients<1,true,false> - (temp1, temp2); - fe_eval.template apply_values<2,true,false> - (temp2, fe_eval.gradients_quad[c][1]); - } + template + inline + void + do_evaluate (FEEval &fe_eval, + const bool evaluate_val, + const bool evaluate_grad, + const bool evaluate_lapl, + internal::int2type<3>) + { + AssertDimension(FEEval::dimension, 3); + Assert (fe_eval.cell != numbers::invalid_unsigned_int, + ExcNotInitialized()); + Assert (fe_eval.dof_values_initialized == true, + internal::ExcAccessToUninitializedField()); - if (evaluate_lapl == true) - { - // grad yz - if (evaluate_grad == false) - fe_eval.template apply_gradients<1,true,false> - (temp1, temp2); - fe_eval.template apply_gradients<2,true,false> - (temp2, fe_eval.hessians_quad[c][5]); + const unsigned int temp_size = FEEval::dofs_per_cell > FEEval::n_q_points ? + FEEval::dofs_per_cell : FEEval::n_q_points; + const unsigned int n_components = FEEval::n_components; - // grad yy - fe_eval.template apply_hessians<1,true,false> - (temp1, temp2); - fe_eval.template apply_values<2,true,false> - (temp2, fe_eval.hessians_quad[c][1]); - } + for (unsigned int c=0; c temp1[temp_size]; + VectorizedArray temp2[temp_size]; - // grad z: can use the values applied in x direction stored in temp1 + if (evaluate_grad == true) + { + // grad x + fe_eval.template apply_gradients<0,true,false> + (fe_eval.values_dofs[c], temp1); fe_eval.template apply_values<1,true,false> - (temp1, temp2); - if (evaluate_grad == true) - fe_eval.template apply_gradients<2,true,false> - (temp2, fe_eval.gradients_quad[c][2]); - - // grad zz: can use the values applied in x and y direction stored - // in temp2 - if (evaluate_lapl == true) - fe_eval.template apply_hessians<2,true,false> - (temp2, fe_eval.hessians_quad[c][2]); - - // val: can use the values applied in x & y direction stored in temp2 - if (evaluate_val == true) - fe_eval.template apply_values<2,true,false> - (temp2, fe_eval.values_quad[c]); - - break; - - case 2: + (temp1, temp2); + fe_eval.template apply_values<2,true,false> + (temp2, fe_eval.gradients_quad[c][0]); + } - // grad x - if (evaluate_grad == true) + if (evaluate_lapl == true) + { + // grad xz + if (evaluate_grad == false) { fe_eval.template apply_gradients<0,true,false> - (fe_eval.values_dofs[c], temp1); - fe_eval.template apply_values<1,true,false> - (temp1, fe_eval.gradients_quad[c][0]); - } - if (evaluate_lapl == true) - { - // grad xy - if (evaluate_grad == false) - fe_eval.template apply_gradients<0,true,false> (fe_eval.values_dofs[c], temp1); - fe_eval.template apply_gradients<1,true,false> - (temp1, fe_eval.hessians_quad[c][2]); - - // grad xx - fe_eval.template apply_hessians<0,true,false> - (fe_eval.values_dofs[c], temp1); fe_eval.template apply_values<1,true,false> - (temp1, fe_eval.hessians_quad[c][0]); + (temp1, temp2); } + fe_eval.template apply_gradients<2,true,false> + (temp2, fe_eval.hessians_quad[c][4]); + + // grad xy + fe_eval.template apply_gradients<1,true,false> + (temp1, temp2); + fe_eval.template apply_values<2,true,false> + (temp2, fe_eval.hessians_quad[c][3]); + + // grad xx + fe_eval.template apply_hessians<0,true,false> + (fe_eval.values_dofs[c], temp1); + fe_eval.template apply_values<1,true,false> + (temp1, temp2); + fe_eval.template apply_values<2,true,false> + (temp2, fe_eval.hessians_quad[c][0]); + } - // grad y - fe_eval.template apply_values<0,true,false> - (fe_eval.values_dofs[c], temp1); - if (evaluate_grad == true) + // grad y + fe_eval.template apply_values<0,true,false> + (fe_eval.values_dofs[c], temp1); + if (evaluate_grad == true) + { + fe_eval.template apply_gradients<1,true,false> + (temp1, temp2); + fe_eval.template apply_values<2,true,false> + (temp2, fe_eval.gradients_quad[c][1]); + } + + if (evaluate_lapl == true) + { + // grad yz + if (evaluate_grad == false) fe_eval.template apply_gradients<1,true,false> - (temp1, fe_eval.gradients_quad[c][1]); + (temp1, temp2); + fe_eval.template apply_gradients<2,true,false> + (temp2, fe_eval.hessians_quad[c][5]); // grad yy - if (evaluate_lapl == true) - fe_eval.template apply_hessians<1,true,false> - (temp1, fe_eval.hessians_quad[c][1]); - - // val: can use values applied in x - if (evaluate_val == true) - fe_eval.template apply_values<1,true,false> - (temp1, fe_eval.values_quad[c]); - - break; - - case 1: - if (evaluate_val == true) - fe_eval.template apply_values<0,true,false> - (fe_eval.values_dofs[c], fe_eval.values_quad[c]); - if (evaluate_grad == true) - fe_eval.template apply_gradients<0,true,false> - (fe_eval.values_dofs[c], fe_eval.gradients_quad[c][0]); - if (evaluate_lapl == true) - fe_eval.template apply_hessians<0,true,false> - (fe_eval.values_dofs[c], fe_eval.hessians_quad[c][0]); - break; - - default: - Assert (false, ExcNotImplemented()); + fe_eval.template apply_hessians<1,true,false> + (temp1, temp2); + fe_eval.template apply_values<2,true,false> + (temp2, fe_eval.hessians_quad[c][1]); } + + // grad z: can use the values applied in x direction stored in temp1 + fe_eval.template apply_values<1,true,false> + (temp1, temp2); + if (evaluate_grad == true) + fe_eval.template apply_gradients<2,true,false> + (temp2, fe_eval.gradients_quad[c][2]); + + // grad zz: can use the values applied in x and y direction stored + // in temp2 + if (evaluate_lapl == true) + fe_eval.template apply_hessians<2,true,false> + (temp2, fe_eval.hessians_quad[c][2]); + + // val: can use the values applied in x & y direction stored in temp2 + if (evaluate_val == true) + fe_eval.template apply_values<2,true,false> + (temp2, fe_eval.values_quad[c]); } #ifdef DEBUG @@ -5059,7 +5125,47 @@ namespace internal void do_integrate (FEEval &fe_eval, const bool integrate_val, - const bool integrate_grad) + const bool integrate_grad, + internal::int2type<1>) + { + Assert (fe_eval.cell != numbers::invalid_unsigned_int, ExcNotInitialized()); + if (integrate_val == true) + Assert (fe_eval.values_quad_submitted == true, + ExcAccessToUninitializedField()); + if (integrate_grad == true) + Assert (fe_eval.gradients_quad_submitted == true, + ExcAccessToUninitializedField()); + + const unsigned int n_components = FEEval::n_components; + + for (unsigned int c=0; c + (fe_eval.gradients_quad[c][0], fe_eval.values_dofs[c]); + if (integrate_val == true) + { + if (integrate_grad == true) + fe_eval.template apply_values<0,false,true> + (fe_eval.values_quad[c], fe_eval.values_dofs[c]); + else + fe_eval.template apply_values<0,false,false> + (fe_eval.values_quad[c], fe_eval.values_dofs[c]); + } + } + +#ifdef DEBUG + fe_eval.dof_values_initialized = true; +#endif + } + + template + inline + void + do_integrate (FEEval &fe_eval, + const bool integrate_val, + const bool integrate_grad, + internal::int2type<2>) { Assert (fe_eval.cell != numbers::invalid_unsigned_int, ExcNotInitialized()); if (integrate_val == true) @@ -5072,109 +5178,108 @@ namespace internal const unsigned int temp_size = FEEval::dofs_per_cell > FEEval::n_q_points ? FEEval::dofs_per_cell : FEEval::n_q_points; const unsigned int n_components = FEEval::n_components; - const unsigned int dim = FEEval::dimension; - for (unsigned int c=0; c temp1[temp_size]; VectorizedArray temp2[temp_size]; - switch (dim) + // val + if (integrate_val == true) + fe_eval.template apply_values<0,false,false> + (fe_eval.values_quad[c], temp1); + if (integrate_grad == true) { - case 3: - + //grad x if (integrate_val == true) - { - // val - fe_eval.template apply_values<0,false,false> - (fe_eval.values_quad[c], temp1); - } - if (integrate_grad == true) - { - // grad x: can sum to temporary value in temp1 - if (integrate_val == true) - fe_eval.template apply_gradients<0,false,true> - (fe_eval.gradients_quad[c][0], temp1); - else - fe_eval.template apply_gradients<0,false,false> - (fe_eval.gradients_quad[c][0], temp1); - } - if (integrate_val || integrate_grad) - fe_eval.template apply_values<1,false,false> - (temp1, temp2); - if (integrate_grad == true) - { - // grad y: can sum to temporary x value in temp2 - fe_eval.template apply_values<0,false,false> - (fe_eval.gradients_quad[c][1], temp1); - fe_eval.template apply_gradients<1,false,true> - (temp1, temp2); - } - if (integrate_val || integrate_grad) - fe_eval.template apply_values<2,false,false> - (temp2, fe_eval.values_dofs[c]); - if (integrate_grad == true) - { - // grad z: can sum to temporary x and y value in output - fe_eval.template apply_values<0,false,false> - (fe_eval.gradients_quad[c][2], temp1); - fe_eval.template apply_values<1,false,false> - (temp1, temp2); - fe_eval.template apply_gradients<2,false,true> - (temp2, fe_eval.values_dofs[c]); - } - - break; + fe_eval.template apply_gradients<0,false,true> + (fe_eval.gradients_quad[c][0], temp1); + else + fe_eval.template apply_gradients<0,false,false> + (fe_eval.gradients_quad[c][0], temp1); + } + if (integrate_val || integrate_grad) + fe_eval.template apply_values<1,false,false> + (temp1, fe_eval.values_dofs[c]); + if (integrate_grad == true) + { + // grad y + fe_eval.template apply_values<0,false,false> + (fe_eval.gradients_quad[c][1], temp1); + fe_eval.template apply_gradients<1,false,true> + (temp1, fe_eval.values_dofs[c]); + } + } - case 2: +#ifdef DEBUG + fe_eval.dof_values_initialized = true; +#endif + } - // val - if (integrate_val == true) - fe_eval.template apply_values<0,false,false> - (fe_eval.values_quad[c], temp1); - if (integrate_grad == true) - { - //grad x - if (integrate_val == true) - fe_eval.template apply_gradients<0,false,true> - (fe_eval.gradients_quad[c][0], temp1); - else - fe_eval.template apply_gradients<0,false,false> - (fe_eval.gradients_quad[c][0], temp1); - } - if (integrate_val || integrate_grad) - fe_eval.template apply_values<1,false,false> - (temp1, fe_eval.values_dofs[c]); - if (integrate_grad == true) - { - // grad y - fe_eval.template apply_values<0,false,false> - (fe_eval.gradients_quad[c][1], temp1); - fe_eval.template apply_gradients<1,false,true> - (temp1, fe_eval.values_dofs[c]); - } + template + inline + void + do_integrate (FEEval &fe_eval, + const bool integrate_val, + const bool integrate_grad, + internal::int2type<3>) + { + Assert (fe_eval.cell != numbers::invalid_unsigned_int, ExcNotInitialized()); + if (integrate_val == true) + Assert (fe_eval.values_quad_submitted == true, + ExcAccessToUninitializedField()); + if (integrate_grad == true) + Assert (fe_eval.gradients_quad_submitted == true, + ExcAccessToUninitializedField()); - break; + const unsigned int temp_size = FEEval::dofs_per_cell > FEEval::n_q_points ? + FEEval::dofs_per_cell : FEEval::n_q_points; + const unsigned int n_components = FEEval::n_components; - case 1: + for (unsigned int c=0; c temp1[temp_size]; + VectorizedArray temp2[temp_size]; - if (integrate_grad == true) - fe_eval.template apply_gradients<0,false,false> - (fe_eval.gradients_quad[c][0], fe_eval.values_dofs[c]); + if (integrate_val == true) + { + // val + fe_eval.template apply_values<0,false,false> + (fe_eval.values_quad[c], temp1); + } + if (integrate_grad == true) + { + // grad x: can sum to temporary value in temp1 if (integrate_val == true) - { - if (integrate_grad == true) - fe_eval.template apply_values<0,false,true> - (fe_eval.values_quad[c], fe_eval.values_dofs[c]); - else - fe_eval.template apply_values<0,false,false> - (fe_eval.values_quad[c], fe_eval.values_dofs[c]); - } - break; - - default: - Assert (false, ExcNotImplemented()); + fe_eval.template apply_gradients<0,false,true> + (fe_eval.gradients_quad[c][0], temp1); + else + fe_eval.template apply_gradients<0,false,false> + (fe_eval.gradients_quad[c][0], temp1); + } + if (integrate_val || integrate_grad) + fe_eval.template apply_values<1,false,false> + (temp1, temp2); + if (integrate_grad == true) + { + // grad y: can sum to temporary x value in temp2 + fe_eval.template apply_values<0,false,false> + (fe_eval.gradients_quad[c][1], temp1); + fe_eval.template apply_gradients<1,false,true> + (temp1, temp2); + } + if (integrate_val || integrate_grad) + fe_eval.template apply_values<2,false,false> + (temp2, fe_eval.values_dofs[c]); + if (integrate_grad == true) + { + // grad z: can sum to temporary x and y value in output + fe_eval.template apply_values<0,false,false> + (fe_eval.gradients_quad[c][2], temp1); + fe_eval.template apply_values<1,false,false> + (temp1, temp2); + fe_eval.template apply_gradients<2,false,true> + (temp2, fe_eval.values_dofs[c]); } } @@ -5196,7 +5301,8 @@ FEEvaluationGeneral const bool evaluate_grad, const bool evaluate_lapl) { - internal::do_evaluate (*this, evaluate_val, evaluate_grad, evaluate_lapl); + internal::do_evaluate (*this, evaluate_val, evaluate_grad, evaluate_lapl, + internal::int2type()); } @@ -5209,7 +5315,8 @@ FEEvaluationGeneral ::integrate (const bool integrate_val, const bool integrate_grad) { - internal::do_integrate (*this, integrate_val, integrate_grad); + internal::do_integrate (*this, integrate_val, integrate_grad, + internal::int2type()); } @@ -5419,7 +5526,8 @@ FEEvaluation const bool evaluate_grad, const bool evaluate_lapl) { - internal::do_evaluate (*this, evaluate_val, evaluate_grad, evaluate_lapl); + internal::do_evaluate (*this, evaluate_val, evaluate_grad, evaluate_lapl, + internal::int2type()); } @@ -5431,7 +5539,8 @@ void FEEvaluation ::integrate (bool integrate_val,bool integrate_grad) { - internal::do_integrate (*this, integrate_val, integrate_grad); + internal::do_integrate (*this, integrate_val, integrate_grad, + internal::int2type()); } -- 2.39.5