From: Shiva Rudraraju Date: Tue, 4 Nov 2014 15:42:47 +0000 (-0500) Subject: Resolve two conflicting partial specializations by adding a third, most specialized... X-Git-Tag: v8.2.0-rc1~89^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=36406765c4c1a3541204af736bcfbad740390c27;p=dealii.git Resolve two conflicting partial specializations by adding a third, most specialized case. --- diff --git a/doc/news/changes.h b/doc/news/changes.h index e7bc34b50b..dde13a8e69 100644 --- a/doc/news/changes.h +++ b/doc/news/changes.h @@ -325,6 +325,13 @@ inconvenience this causes.

Specific improvements

    +
  1. Fixed: Using the FEEvaluation framework did not work for + scalar elements in 1d because there were conflicting partial + specializations. This is now fixed. +
    + (Shiva Rudraraju, 2014/11/04) +
    2. +
  2. New: There is now a macro DEAL_II_VERSION_GTE that can be used to test whether the deal.II version is greater than or equal a particular version number. This is useful if you diff --git a/include/deal.II/matrix_free/fe_evaluation.h b/include/deal.II/matrix_free/fe_evaluation.h index bd8ed69797..c143e2f5ad 100644 --- a/include/deal.II/matrix_free/fe_evaluation.h +++ b/include/deal.II/matrix_free/fe_evaluation.h @@ -1224,6 +1224,162 @@ protected: }; +/** + * This class provides access to the data fields of the FEEvaluation + * classes. Partial specialization for scalar fields in 1d that defines access with + * simple data fields, i.e., scalars for the values and Tensor<1,1> for the + * gradients. + * + * @author Katharina Kormann and Martin Kronbichler, 2010, 2011, Shiva Rudraraju, 2014 + */ +template +class FEEvaluationAccess<1,1,Number> : public FEEvaluationBase<1,1,Number> +{ +public: + static const int dim = 1; + typedef Number number_type; + typedef VectorizedArray value_type; + typedef Tensor<1,dim,VectorizedArray > gradient_type; + static const unsigned int dimension = dim; + typedef FEEvaluationBase BaseClass; + + /** + * Returns 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 + * function of the given index. + */ + value_type get_dof_value (const unsigned int dof) const; + + /** + * Write a value to the field containing the degrees of freedom with + * component @p dof. Access to the same field as through @p get_dof_value. + */ + void submit_dof_value (const value_type val_in, + const unsigned int dof); + + /** + * Returns the value of a finite element function at quadrature point number + * @p q_point after a call to @p evaluate(true,...), or the value that has + * been stored there with a call to @p submit_value. 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. + */ + value_type get_value (const unsigned int q_point) const; + + /** + * Write a value to the field containing the values on quadrature points + * with component @p q_point. Access to the same field as through @p + * get_value. If applied before the function @p integrate(true,...) is + * called, this specifies the value which is tested by all basis function on + * the current cell and integrated over. + */ + void submit_value (const value_type val_in, + const unsigned int q_point); + + /** + * Returns the gradient of a finite element function at quadrature point + * number @p q_point after a call to @p evaluate(...,true,...), or the value + * that has been stored there with a call to @p submit_gradient. + */ + gradient_type get_gradient (const unsigned int q_point) const; + + /** + * Write a contribution that is tested by the gradient to the field + * containing the values on quadrature points with component @p + * q_point. Access to the same field as through @p get_gradient. If applied + * before the function @p integrate(...,true) is called, this specifies what + * is tested by all basis function gradients on the current cell and + * integrated over. + */ + void submit_gradient(const gradient_type grad_in, + const unsigned int q_point); + + /** + * Returns the Hessian of a finite element function at quadrature point + * number @p q_point after a call to @p evaluate(...,true). If only the + * diagonal part of the Hessian or its trace, the Laplacian, are needed, use + * the respective functions below. + */ + Tensor<2,1,VectorizedArray > + get_hessian (unsigned int q_point) const; + + /** + * Returns the diagonal of the Hessian of a finite element function at + * quadrature point number @p q_point after a call to @p evaluate(...,true). + */ + gradient_type get_hessian_diagonal (const unsigned int q_point) const; + + /** + * Returns the Laplacian of a finite element function at quadrature point + * number @p q_point after a call to @p evaluate(...,true). + */ + value_type get_laplacian (const unsigned int q_point) const; + + /** + * Takes values on quadrature points, multiplies by the Jacobian determinant + * and quadrature weights (JxW) and sums the values for all quadrature + * points on the cell. The result is a scalar, representing the integral + * over the function over the cell. If a vector-element is used, the + * resulting components are still separated. Moreover, if vectorization is + * enabled, the integral values of several cells are represented together. + */ + value_type integrate_value () const; + +protected: + /** + * Constructor. Made protected to avoid initialization in user code. Takes + * all data stored in MatrixFree. If applied to problems with more than one + * finite element or more than one quadrature formula selected during + * construction of @p matrix_free, @p fe_no and @p quad_no allow to select + * the appropriate components. + */ + FEEvaluationAccess (const MatrixFree<1,Number> &matrix_free, + const unsigned int fe_no, + const unsigned int quad_no, + const unsigned int fe_degree, + const unsigned int n_q_points); + + /** + * Constructor that comes with reduced functionality and works similar as + * FEValues. The user has to provide a structure of type MappingFEEvaluation + * and a DoFHandler in order to allow for reading out the finite element + * data. It uses the data provided by dof_handler.get_fe(). If the element + * is vector-valued, the optional argument allows to specify the index of + * the base element (as long as the element is primitive, non-primitive are + * not supported currently). + * + * With initialization from a FEValues object, no call to a reinit method of + * this class is necessary. Instead, it is enough if the geometry is + * initialized to a given cell iterator. It can also read from or write to + * vectors in the standard way for DoFHandler::active_cell_iterator + * types (which is less efficient with MPI since index translation has to be + * done), but of course only for one cell at a time. Hence, a kernel using + * this method does not vectorize over several elements (which is most + * efficient for vector operations), but only possibly within the element if + * the evaluate/integrate routines are combined (e.g. for computing cell + * matrices). + * With this initialization, no call to a reinit method of this + * class. Instead, it is enough if the geometry is initialized to a given + * cell iterator. Moreover, beware that a kernel using this method does not + * vectorize over several elements (which is most efficient for vector + * operations), but only possibly within the element if the + * evaluate/integrate routines are combined (e.g. for matrix assembly). + */ + FEEvaluationAccess (const MappingFEEvaluation<1,Number> &geometry, + const DoFHandler<1> &dof_handler, + const unsigned int first_selected_component = 0); + + /** + * Copy constructor + */ + FEEvaluationAccess (const FEEvaluationAccess &other); +}; + + /** * The class that provides all functions necessary to evaluate functions at @@ -4076,6 +4232,230 @@ FEEvaluationAccess } +/*-------------------- FEEvaluationAccess scalar for 1d ----------------------------*/ + + +template +inline +FEEvaluationAccess<1,1,Number> +::FEEvaluationAccess (const MatrixFree<1,Number> &data_in, + const unsigned int fe_no, + const unsigned int quad_no_in, + const unsigned int fe_degree, + const unsigned int n_q_points) + : + FEEvaluationBase <1,1,Number> + (data_in, fe_no, quad_no_in, fe_degree, n_q_points) +{} + + + +template +inline +FEEvaluationAccess<1,1,Number> +::FEEvaluationAccess (const MappingFEEvaluation<1,Number> &geometry, + const DoFHandler<1> &dof_handler, + const unsigned int first_selected_component) + : + FEEvaluationBase <1,1,Number> (geometry, dof_handler, first_selected_component) +{} + + + +template +inline +FEEvaluationAccess<1,1,Number> +::FEEvaluationAccess (const FEEvaluationAccess &other) + : + FEEvaluationBase <1,1,Number>(other) +{} + + + +template +inline +VectorizedArray +FEEvaluationAccess<1,1,Number> +::get_dof_value (const unsigned int dof) const +{ + AssertIndexRange (dof, this->data->dofs_per_cell); + return this->values_dofs[0][dof]; +} + + + +template +inline +VectorizedArray +FEEvaluationAccess<1,1,Number> +::get_value (const unsigned int q_point) const +{ + Assert (this->values_quad_initialized==true, + internal::ExcAccessToUninitializedField()); + AssertIndexRange (q_point, this->data->n_q_points); + return this->values_quad[0][q_point]; +} + + + +template +inline +Tensor<1,1,VectorizedArray > +FEEvaluationAccess<1,1,Number> +::get_gradient (const unsigned int q_point) const +{ + // could use the base class gradient, but that involves too many inefficient + // initialization operations on tensors + + Assert (this->gradients_quad_initialized==true, + internal::ExcAccessToUninitializedField()); + AssertIndexRange (q_point, this->data->n_q_points); + + Tensor<1,1,VectorizedArray > grad_out (false); + + // Cartesian cell + if (this->cell_type == internal::MatrixFreeFunctions::cartesian) + { + grad_out[0] = (this->gradients_quad[0][0][q_point] * + this->cartesian_data[0][0]); + } + // cell with general/constant Jacobian + else + { + const Tensor<2,1,VectorizedArray > &jac = + this->cell_type == internal::MatrixFreeFunctions::general ? + this->jacobian[q_point] : this->jacobian[0]; + + grad_out[0] = (jac[0][0] * this->gradients_quad[0][0][q_point]); + } + return grad_out; +} + + + +template +inline +Tensor<2,1,VectorizedArray > +FEEvaluationAccess<1,1,Number> +::get_hessian (const unsigned int q_point) const +{ + return BaseClass::get_hessian(q_point)[0]; +} + + + +template +inline +Tensor<1,1,VectorizedArray > +FEEvaluationAccess<1,1,Number> +::get_hessian_diagonal (const unsigned int q_point) const +{ + return BaseClass::get_hessian_diagonal(q_point)[0]; +} + + + +template +inline +VectorizedArray +FEEvaluationAccess<1,1,Number> +::get_laplacian (const unsigned int q_point) const +{ + return BaseClass::get_laplacian(q_point)[0]; +} + + + +template +inline +void +FEEvaluationAccess<1,1,Number> +::submit_dof_value (const VectorizedArray val_in, + const unsigned int dof) +{ +#ifdef DEBUG + this->dof_values_initialized = true; + AssertIndexRange (dof, this->data->dofs_per_cell); +#endif + this->values_dofs[0][dof] = val_in; +} + + + +template +inline +void +FEEvaluationAccess<1,1,Number> +::submit_value (const VectorizedArray val_in, + const unsigned int q_point) +{ +#ifdef DEBUG + Assert (this->cell != numbers::invalid_unsigned_int, ExcNotInitialized()); + AssertIndexRange (q_point, this->data->n_q_points); + this->values_quad_submitted = true; +#endif + if (this->cell_type == internal::MatrixFreeFunctions::general) + { + const VectorizedArray JxW = this->J_value[q_point]; + this->values_quad[0][q_point] = val_in * JxW; + } + else //if (this->cell_type < internal::MatrixFreeFunctions::general) + { + const VectorizedArray JxW = this->J_value[0] * this->quadrature_weights[q_point]; + this->values_quad[0][q_point] = val_in * JxW; + } +} + + + +template +inline +void +FEEvaluationAccess<1,1,Number> +::submit_gradient (const Tensor<1,1,VectorizedArray > grad_in, + const unsigned int q_point) +{ +#ifdef DEBUG + Assert (this->cell != numbers::invalid_unsigned_int, ExcNotInitialized()); + AssertIndexRange (q_point, this->data->n_q_points); + this->gradients_quad_submitted = true; +#endif + if (this->cell_type == internal::MatrixFreeFunctions::cartesian) + { + const VectorizedArray JxW = this->J_value[0] * this->quadrature_weights[q_point]; + this->gradients_quad[0][0][q_point] = (grad_in[0] * + this->cartesian_data[0][0] * + JxW); + } + // general/affine cell type + else + { + const Tensor<2,1,VectorizedArray > &jac = + this->cell_type == internal::MatrixFreeFunctions::general ? + this->jacobian[q_point] : this->jacobian[0]; + const VectorizedArray JxW = + this->cell_type == internal::MatrixFreeFunctions::general ? + this->J_value[q_point] : this->J_value[0] * this->quadrature_weights[q_point]; + + VectorizedArray new_val = jac[0][0] * grad_in[0]; + for (unsigned int e=1; egradients_quad[0][0][q_point] = new_val * JxW; + } +} + + + +template +inline +VectorizedArray +FEEvaluationAccess<1,1,Number> +::integrate_value () const +{ + return BaseClass::integrate_value()[0]; +} + + namespace internal { diff --git a/tests/matrix_free/fe_evaluation_access_1d.cc b/tests/matrix_free/fe_evaluation_access_1d.cc new file mode 100644 index 0000000000..c4d5267dee --- /dev/null +++ b/tests/matrix_free/fe_evaluation_access_1d.cc @@ -0,0 +1,34 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2014 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + + + +// FEEvaluationAccess<1,1,double> didn't compile because we had +// conflicting partial specializations of this class + +#include "../tests.h" +#include + + +int main () +{ + std::ofstream logfile("output"); + deallog.attach(logfile); + deallog.depth_console(0); + deallog.threshold_double(1.e-10); + + FEEvaluationAccess<1,1,double> *test; // didn't compile before + deallog << "OK" << std::endl; +} diff --git a/tests/matrix_free/fe_evaluation_access_1d.output b/tests/matrix_free/fe_evaluation_access_1d.output new file mode 100644 index 0000000000..0fd8fc12f0 --- /dev/null +++ b/tests/matrix_free/fe_evaluation_access_1d.output @@ -0,0 +1,2 @@ + +DEAL::OK