From: Jiaqi Zhang Date: Mon, 9 Aug 2021 19:47:40 +0000 (+0000) Subject: add get_function_* X-Git-Tag: v9.4.0-rc1~1090^2~1 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=c3fb164c7da3570599c73697b1ade11b4c6793bd;p=dealii.git add get_function_* --- diff --git a/include/deal.II/fe/fe_interface_values.h b/include/deal.II/fe/fe_interface_values.h index 89503d44f6..b3ecab1cc3 100644 --- a/include/deal.II/fe/fe_interface_values.h +++ b/include/deal.II/fe/fe_interface_values.h @@ -55,6 +55,15 @@ namespace FEInterfaceViews * Store a pointer to the FEInterfaceValues instance. */ const FEInterfaceValues *fe_interface; + + /** + * Get the local degree-of-freedom values associated with the internally + * initialized cell interface. + */ + template + void + get_local_dof_values(const InputVector &dof_values, + OutputVector & local_dof_values) const; }; @@ -92,6 +101,45 @@ namespace FEInterfaceViews using third_derivative_type = typename FEValuesViews::Scalar::third_derivative_type; + /** + * An alias for the data type of the product of a @p Number and the + * values of the view this class provides. This is the data type of + * scalar components of a finite element field whose degrees of + * freedom are described by a vector with elements of type @p Number. + */ + template + using solution_value_type = typename ProductType::type; + + /** + * An alias for the data type of the product of a @p Number and the + * gradients of the view this class provides. This is the data type of + * scalar components of a finite element field whose degrees of + * freedom are described by a vector with elements of type @p Number. + */ + template + using solution_gradient_type = + typename ProductType::type; + + /** + * An alias for the data type of the product of a @p Number and the + * Hessians of the view this class provides. This is the data type of + * scalar components of a finite element field whose degrees of + * freedom are described by a vector with elements of type @p Number. + */ + template + using solution_hessian_type = + typename ProductType::type; + + /** + * An alias for the data type of the product of a @p Number and the + * third derivatives of the view this class provides. This is the data type + * of scalar components of a finite element field whose degrees of + * freedom are described by a vector with elements of type @p Number. + */ + template + using solution_third_derivative_type = + typename ProductType::type; + /** * Constructor for an object that represents a single scalar component */ @@ -273,6 +321,304 @@ namespace FEInterfaceViews jump_3rd_derivative(const unsigned int interface_dof_index, const unsigned int q_point) const; + /** + * Return the values of the selected scalar component of the finite + * element function characterized by fe_function at the + * quadrature points of the cell interface selected the last time + * the reinit function of the FEInterfaceValues object was called. + * + * The argument @p here_or_there selects between the value on cell 0 (here, @p true) + * and cell 1 (there, @p false). You can also interpret it as "upstream" (@p true) + * and "downstream" (@p false) as defined by the direction of the normal + * vector in this quadrature point. If @p here_or_there is true, the values + * from the first cell of the interface is used. + * + * The data type stored by the output vector must be what you get when you + * multiply the values of shape functions (i.e., @p value_type) times the + * type used to store the values of the unknowns $U_j$ of your finite + * element vector $U$ (represented by the @p fe_function argument). + * + * @dealiiRequiresUpdateFlags{update_values} + */ + template + void + get_function_values( + const bool here_or_there, + const InputVector &fe_function, + std::vector> + &values) const; + + /** + * Same as above, but using a vector of local degree-of-freedom values. In + * other words, instead of extracting the nodal values of the degrees of + * freedom located on the current cell interface from a global vector + * associated with a DoFHandler object (as the function above does), this + * function instead takes these local nodal values through its first + * argument. + * + * @param[in] here_or_there Same as the one in the above function. + * + * @param[in] local_dof_values A vector of local nodal values. This vector + * must have a length equal to number of DoFs on the current cell, and + * must be ordered in the same order as degrees of freedom are numbered on + * the reference cell. + * + * @param[out] values A vector of values of the given finite element field, + * at the quadrature points on the current object. + * + * @tparam InputVector The @p InputVector type must allow creation + * of an ArrayView object from it; this is satisfied by the + * `std::vector` class, among others. + */ + template + void + get_function_values_from_local_dof_values( + const bool here_or_there, + const InputVector &local_dof_values, + std::vector> + &values) const; + + /** + * Return the jump in the values of the selected scalar component of the + * finite element function characterized by fe_function at the + * quadrature points of the cell interface selected the last time + * the reinit function of the FEInterfaceValues object was called. + * + * The data type stored by the output vector must be what you get when you + * multiply the values of shape functions (i.e., @p value_type) times the + * type used to store the values of the unknowns $U_j$ of your finite + * element vector $U$ (represented by the @p fe_function argument). + * + * @dealiiRequiresUpdateFlags{update_values} + */ + template + void + get_jump_in_function_values( + const InputVector &fe_function, + std::vector> + &values) const; + + /** + * This function relates to get_jump_in_function_values() in the same way + * as get_function_values_from_local_dof_values() relates to + * get_function_values(). See the documentation of + * get_function_values_from_local_dof_values() for more + * information. + */ + template + void + get_jump_in_function_values_from_local_dof_values( + const InputVector &local_dof_values, + std::vector> + &values) const; + + /** + * Return the average of the values of the selected scalar component of the + * finite element function characterized by fe_function at the + * quadrature points of the cell interface selected the last time + * the reinit function of the FEInterfaceValues object was called. + * + * The data type stored by the output vector must be what you get when you + * multiply the values of shape functions (i.e., @p value_type) times the + * type used to store the values of the unknowns $U_j$ of your finite + * element vector $U$ (represented by the @p fe_function argument). + * + * @dealiiRequiresUpdateFlags{update_values} + */ + template + void + get_average_of_function_values( + const InputVector &fe_function, + std::vector> + &values) const; + + /** + * This function relates to get_average_of_function_values() in the same way + * as get_function_values_from_local_dof_values() relates to + * get_function_values(). See the documentation of + * get_function_values_from_local_dof_values() for more + * information. + */ + template + void + get_average_of_function_values_from_local_dof_values( + const InputVector &local_dof_values, + std::vector> + &values) const; + + /** + * Return the jump in the gradients of the selected scalar components of the + * finite element function characterized by fe_function at the + * quadrature points of the cell interface selected the last time + * the reinit function of the FEInterfaceValues object was called. + * + * The data type stored by the output vector must be what you get when you + * multiply the gradients of shape functions (i.e., @p gradient_type) + * times the type used to store the values of the unknowns $U_j$ of your + * finite element vector $U$ (represented by the @p fe_function argument). + * + * @dealiiRequiresUpdateFlags{update_gradients} + */ + template + void + get_jump_in_function_gradients( + const InputVector &fe_function, + std::vector> + &gradients) const; + + /** + * This function relates to get_jump_in_function_gradients() in the same way + * as get_function_values_from_local_dof_values() relates to + * get_function_values(). See the documentation of + * get_function_values_from_local_dof_values() for more + * information. + */ + template + void + get_jump_in_function_gradients_from_local_dof_values( + const InputVector &local_dof_values, + std::vector> + &gradients) const; + + /** + * Return the average of the gradients of the selected scalar components of + * the finite element function characterized by fe_function at the + * quadrature points of the cell interface selected the last time + * the reinit function of the FEInterfaceValues object was called. + * + * The data type stored by the output vector must be what you get when you + * multiply the gradients of shape functions (i.e., @p gradient_type) + * times the type used to store the values of the unknowns $U_j$ of your + * finite element vector $U$ (represented by the @p fe_function argument). + * + * @dealiiRequiresUpdateFlags{update_gradients} + */ + template + void + get_average_of_function_gradients( + const InputVector &fe_function, + std::vector> + &gradients) const; + + /** + * This function relates to get_average_of_function_gradients() in the same + * way as get_function_values_from_local_dof_values() relates to + * get_function_values(). See the documentation of + * get_function_values_from_local_dof_values() for more + * information. + */ + template + void + get_average_of_function_gradients_from_local_dof_values( + const InputVector &local_dof_values, + std::vector> + &gradients) const; + + /** + * Return the average of the Hessians of the selected scalar component of + * the finite element function characterized by fe_function at the + * quadrature points of the cell, face or subface selected the last time + * the reinit function of the FEInterfaceValues object was called. + * + * The data type stored by the output vector must be what you get when you + * multiply the Hessians of shape functions (i.e., @p hessian_type) times + * the type used to store the values of the unknowns $U_j$ of your finite + * element vector $U$ (represented by the @p fe_function argument). + * + * @dealiiRequiresUpdateFlags{update_hessians} + */ + template + void + get_average_of_function_hessians( + const InputVector &fe_function, + std::vector> + &hessians) const; + + /** + * This function relates to get_average_of_function_hessians() in the same + * way as get_function_values_from_local_dof_values() relates to + * get_function_values(). See the documentation of + * get_function_values_from_local_dof_values() for more + * information. + */ + template + void + get_average_of_function_hessians_from_local_dof_values( + const InputVector &local_dof_values, + std::vector> + &hessians) const; + + /** + * Return the jump in the Hessians of the selected scalar component of the + * finite element function characterized by fe_function at the + * quadrature points of the cell, face or subface selected the last time + * the reinit function of the FEInterfaceValues object was called. + * + * The data type stored by the output vector must be what you get when you + * multiply the Hessians of shape functions (i.e., @p hessian_type) times + * the type used to store the values of the unknowns $U_j$ of your finite + * element vector $U$ (represented by the @p fe_function argument). + * + * @dealiiRequiresUpdateFlags{update_hessians} + */ + template + void + get_jump_in_function_hessians( + const InputVector &fe_function, + std::vector> + &hessians) const; + + /** + * This function relates to get_jump_in_function_hessians() in the same way + * as get_function_values_from_local_dof_values() relates to + * get_function_values(). See the documentation of + * get_function_values_from_local_dof_values() for more + * information. + */ + template + void + get_jump_in_function_hessians_from_local_dof_values( + const InputVector &local_dof_values, + std::vector> + &hessians) const; + + /** + * Return the jump in the third derivatives of the selected scalar component + * of the finite element function characterized by fe_function at + * the quadrature points of the cell, face or subface selected the last time + * the reinit function of the FEInterfaceValues object was called. + * + * The data type stored by the output vector must be what you get when you + * multiply the third derivatives of shape functions (i.e., @p + * third_derivative_type) times the type used to store the values of the + * unknowns $U_j$ of your finite element vector $U$ (represented by the @p + * fe_function argument). + * + * @dealiiRequiresUpdateFlags{update_third_derivatives} + */ + template + void + get_jump_in_function_third_derivatives( + const InputVector &fe_function, + std::vector< + solution_third_derivative_type> + &third_derivatives) const; + + /** + * This function relates to get_jump_in_function_third_derivatives() in the + * same way as get_function_values_from_local_dof_values() relates + * to get_function_values(). See the documentation of + * get_function_values_from_local_dof_values() for more + * information. + */ + template + void + get_jump_in_function_third_derivatives_from_local_dof_values( + const InputVector &local_dof_values, + std::vector< + solution_third_derivative_type> + &third_derivatives) const; + private: /** * The extractor for this view. @@ -482,6 +828,344 @@ namespace FEInterfaceViews jump_3rd_derivative(const unsigned int interface_dof_index, const unsigned int q_point) const; + /** + * An alias for the data type of the product of a @p Number and the + * values of the view this class provides. This is the data type of + * vector components of a finite element field whose degrees of + * freedom are described by a vector with elements of type @p Number. + */ + template + using solution_value_type = typename ProductType::type; + + /** + * An alias for the data type of the product of a @p Number and the + * gradients of the view this class provides. This is the data type of + * vector components of a finite element field whose degrees of + * freedom are described by a vector with elements of type @p Number. + */ + template + using solution_gradient_type = + typename ProductType::type; + + /** + * An alias for the data type of the product of a @p Number and the + * Hessians of the view this class provides. This is the data type of + * vector components of a finite element field whose degrees of + * freedom are described by a vector with elements of type @p Number. + */ + template + using solution_hessian_type = + typename ProductType::type; + + /** + * An alias for the data type of the product of a @p Number and the + * third derivatives of the view this class provides. This is the data type + * of vector components of a finite element field whose degrees of + * freedom are described by a vector with elements of type @p Number. + */ + template + using solution_third_derivative_type = + typename ProductType::type; + + /** + * Return the values of the selected vector component of the finite + * element function characterized by fe_function at the + * quadrature points of the cell interface selected the last time + * the reinit function of the FEInterfaceValues object was called. + * + * The argument @p here_or_there selects between the value on cell 0 (here, @p true) + * and cell 1 (there, @p false). You can also interpret it as "upstream" (@p true) + * and "downstream" (@p false) as defined by the direction of the normal + * vector in this quadrature point. If @p here_or_there is true, the values + * from the first cell of the interface is used. + * + * The data type stored by the output vector must be what you get when you + * multiply the values of shape functions (i.e., @p value_type) times the + * type used to store the values of the unknowns $U_j$ of your finite + * element vector $U$ (represented by the @p fe_function argument). + * + * @dealiiRequiresUpdateFlags{update_values} + */ + template + void + get_function_values( + const bool here_or_there, + const InputVector &fe_function, + std::vector> + &values) const; + + /** + * Same as above, but using a vector of local degree-of-freedom values. In + * other words, instead of extracting the nodal values of the degrees of + * freedom located on the current cell interface from a global vector + * associated with a DoFHandler object (as the function above does), this + * function instead takes these local nodal values through its first + * argument. + * + * @param[in] here_or_there Same as the one in the above function. + * + * @param[in] local_dof_values A vector of local nodal values. This vector + * must have a length equal to number of DoFs on the current cell, and + * must be ordered in the same order as degrees of freedom are numbered on + * the reference cell. + * + * @param[out] values A vector of values of the given finite element field, + * at the quadrature points on the current object. + * + * @tparam InputVector The @p InputVector type must allow creation + * of an ArrayView object from it; this is satisfied by the + * `std::vector` class, among others. + */ + template + void + get_function_values_from_local_dof_values( + const bool here_or_there, + const InputVector &local_dof_values, + std::vector> + &values) const; + + /** + * Return the jump in the values of the selected vector component of the + * finite element function characterized by fe_function at the + * quadrature points of the cell interface selected the last time + * the reinit function of the FEInterfaceValues object was called. + * + * The data type stored by the output vector must be what you get when you + * multiply the values of shape functions (i.e., @p value_type) times the + * type used to store the values of the unknowns $U_j$ of your finite + * element vector $U$ (represented by the @p fe_function argument). + * + * @dealiiRequiresUpdateFlags{update_values} + */ + template + void + get_jump_in_function_values( + const InputVector &fe_function, + std::vector> + &values) const; + + /** + * This function relates to get_jump_in_function_values() in the same way + * as get_function_values_from_local_dof_values() relates to + * get_function_values(). See the documentation of + * get_function_values_from_local_dof_values() for more + * information. + */ + template + void + get_jump_in_function_values_from_local_dof_values( + const InputVector &local_dof_values, + std::vector> + &values) const; + + /** + * Return the average of the values of the selected vector component of the + * finite element function characterized by fe_function at the + * quadrature points of the cell interface selected the last time + * the reinit function of the FEInterfaceValues object was called. + * + * The data type stored by the output vector must be what you get when you + * multiply the values of shape functions (i.e., @p value_type) times the + * type used to store the values of the unknowns $U_j$ of your finite + * element vector $U$ (represented by the @p fe_function argument). + * + * @dealiiRequiresUpdateFlags{update_values} + */ + template + void + get_average_of_function_values( + const InputVector &fe_function, + std::vector> + &values) const; + + /** + * This function relates to get_average_of_function_values() in the same way + * as get_function_values_from_local_dof_values() relates to + * get_function_values(). See the documentation of + * get_function_values_from_local_dof_values() for more + * information. + */ + template + void + get_average_of_function_values_from_local_dof_values( + const InputVector &local_dof_values, + std::vector> + &values) const; + + /** + * Return the jump in the gradients of the selected vector components of the + * finite element function characterized by fe_function at the + * quadrature points of the cell interface selected the last time + * the reinit function of the FEInterfaceValues object was called. + * + * The data type stored by the output vector must be what you get when you + * multiply the gradients of shape functions (i.e., @p gradient_type) + * times the type used to store the values of the unknowns $U_j$ of your + * finite element vector $U$ (represented by the @p fe_function argument). + * + * @dealiiRequiresUpdateFlags{update_gradients} + */ + template + void + get_jump_in_function_gradients( + const InputVector &fe_function, + std::vector> + &gradients) const; + + /** + * This function relates to get_jump_in_function_gradients() in the same way + * as get_function_values_from_local_dof_values() relates to + * get_function_values(). See the documentation of + * get_function_values_from_local_dof_values() for more + * information. + */ + template + void + get_jump_in_function_gradients_from_local_dof_values( + const InputVector &local_dof_values, + std::vector> + &gradients) const; + + /** + * Return the average of the gradients of the selected vector components of + * the finite element function characterized by fe_function at the + * quadrature points of the cell interface selected the last time + * the reinit function of the FEInterfaceValues object was called. + * + * The data type stored by the output vector must be what you get when you + * multiply the gradients of shape functions (i.e., @p gradient_type) + * times the type used to store the values of the unknowns $U_j$ of your + * finite element vector $U$ (represented by the @p fe_function argument). + * + * @dealiiRequiresUpdateFlags{update_gradients} + */ + template + void + get_average_of_function_gradients( + const InputVector &fe_function, + std::vector> + &gradients) const; + + /** + * This function relates to get_average_of_function_gradients() in the same + * way as get_function_values_from_local_dof_values() relates to + * get_function_values(). See the documentation of + * get_function_values_from_local_dof_values() for more + * information. + */ + template + void + get_average_of_function_gradients_from_local_dof_values( + const InputVector &local_dof_values, + std::vector> + &gradients) const; + + /** + * Return the average of the Hessians of the selected vector component of + * the finite element function characterized by fe_function at the + * quadrature points of the cell, face or subface selected the last time + * the reinit function of the FEInterfaceValues object was called. + * + * The data type stored by the output vector must be what you get when you + * multiply the Hessians of shape functions (i.e., @p hessian_type) times + * the type used to store the values of the unknowns $U_j$ of your finite + * element vector $U$ (represented by the @p fe_function argument). + * + * @dealiiRequiresUpdateFlags{update_hessians} + */ + template + void + get_average_of_function_hessians( + const InputVector &fe_function, + std::vector> + &hessians) const; + + /** + * This function relates to get_average_of_function_hessians() in the same + * way as get_function_values_from_local_dof_values() relates to + * get_function_values(). See the documentation of + * get_function_values_from_local_dof_values() for more + * information. + */ + template + void + get_average_of_function_hessians_from_local_dof_values( + const InputVector &local_dof_values, + std::vector> + &hessians) const; + + /** + * Return the jump in the Hessians of the selected vector component of the + * finite element function characterized by fe_function at the + * quadrature points of the cell, face or subface selected the last time + * the reinit function of the FEInterfaceValues object was called. + * + * The data type stored by the output vector must be what you get when you + * multiply the Hessians of shape functions (i.e., @p hessian_type) times + * the type used to store the values of the unknowns $U_j$ of your finite + * element vector $U$ (represented by the @p fe_function argument). + * + * @dealiiRequiresUpdateFlags{update_hessians} + */ + template + void + get_jump_in_function_hessians( + const InputVector &fe_function, + std::vector> + &hessians) const; + + /** + * This function relates to get_jump_in_function_hessians() in the same way + * as get_function_values_from_local_dof_values() relates to + * get_function_values(). See the documentation of + * get_function_values_from_local_dof_values() for more + * information. + */ + template + void + get_jump_in_function_hessians_from_local_dof_values( + const InputVector &local_dof_values, + std::vector> + &hessians) const; + + /** + * Return the jump in the third derivatives of the selected vector component + * of the finite element function characterized by fe_function at + * the quadrature points of the cell, face or subface selected the last time + * the reinit function of the FEInterfaceValues object was called. + * + * The data type stored by the output vector must be what you get when you + * multiply the third derivatives of shape functions (i.e., @p + * third_derivative_type) times the type used to store the values of the + * unknowns $U_j$ of your finite element vector $U$ (represented by the @p + * fe_function argument). + * + * @dealiiRequiresUpdateFlags{update_third_derivatives} + */ + template + void + get_jump_in_function_third_derivatives( + const InputVector &fe_function, + std::vector< + solution_third_derivative_type> + &third_derivatives) const; + + /** + * This function relates to get_jump_in_function_third_derivatives() in the + * same way as get_function_values_from_local_dof_values() relates + * to get_function_values(). See the documentation of + * get_function_values_from_local_dof_values() for more + * information. + */ + template + void + get_jump_in_function_third_derivatives_from_local_dof_values( + const InputVector &local_dof_values, + std::vector< + solution_third_derivative_type> + &third_derivatives) const; + + private: /** * The extractor for this view. @@ -1835,6 +2519,24 @@ namespace FEInterfaceViews + template + template + void + Base::get_local_dof_values( + const InputVector &dof_values, + OutputVector & local_dof_values) const + { + const auto &interface_dof_indices = + this->fe_interface->get_interface_dof_indices(); + + AssertDimension(interface_dof_indices.size(), local_dof_values.size()); + + for (unsigned int i = 0; i < interface_dof_indices.size(); ++i) + local_dof_values[i] = dof_values(interface_dof_indices[i]); + } + + + template typename Scalar::value_type Scalar::value(const bool here_or_there, @@ -2125,6 +2827,375 @@ namespace FEInterfaceViews + template + template + void + Scalar::get_function_values_from_local_dof_values( + const bool here_or_there, + const InputVector &local_dof_values, + std::vector> &values) + const + { + Assert(values.size() == this->fe_interface->n_quadrature_points, + ExcNotImplemented()); + + const auto &interface_dof_indices = + this->fe_interface->get_interface_dof_indices(); + + for (const auto dof_index : this->fe_interface->dof_indices()) + for (const auto q_index : this->fe_interface->quadrature_point_indices()) + { + if (dof_index == 0) + values[q_index] = 0.; + + values[q_index] += local_dof_values[dof_index] * + value(here_or_there, dof_index, q_index); + } + } + + + + template + template + void + Scalar::get_function_values( + const bool here_or_there, + const InputVector &fe_function, + std::vector> &values) + const + { + std::vector local_dof_values( + this->fe_interface->n_current_interface_dofs()); + this->get_local_dof_values(fe_function, local_dof_values); + + get_function_values_from_local_dof_values(here_or_there, + local_dof_values, + values); + } + + + + template + template + void + Scalar::get_jump_in_function_values_from_local_dof_values( + const InputVector &local_dof_values, + std::vector> &values) + const + { + Assert(values.size() == this->fe_interface->n_quadrature_points, + ExcNotImplemented()); + + const auto &interface_dof_indices = + this->fe_interface->get_interface_dof_indices(); + + for (const auto dof_index : this->fe_interface->dof_indices()) + for (const auto q_index : this->fe_interface->quadrature_point_indices()) + { + if (dof_index == 0) + values[q_index] = 0.; + + values[q_index] += + local_dof_values[dof_index] * jump_in_values(dof_index, q_index); + } + } + + + + template + template + void + Scalar::get_jump_in_function_values( + const InputVector &fe_function, + std::vector> &values) + const + { + std::vector local_dof_values( + this->fe_interface->n_current_interface_dofs()); + this->get_local_dof_values(fe_function, local_dof_values); + + get_jump_in_function_values_from_local_dof_values(local_dof_values, values); + } + + + + template + template + void + Scalar::get_jump_in_function_gradients_from_local_dof_values( + const InputVector &local_dof_values, + std::vector> + &gradients) const + { + Assert(gradients.size() == this->fe_interface->n_quadrature_points, + ExcNotImplemented()); + + const auto &interface_dof_indices = + this->fe_interface->get_interface_dof_indices(); + + for (const auto dof_index : this->fe_interface->dof_indices()) + for (const auto q_index : this->fe_interface->quadrature_point_indices()) + { + if (dof_index == 0) + gradients[q_index] = 0.; + + gradients[q_index] += + local_dof_values[dof_index] * jump_in_gradients(dof_index, q_index); + } + } + + + + template + template + void + Scalar::get_jump_in_function_gradients( + const InputVector &fe_function, + std::vector> + &gradients) const + { + std::vector local_dof_values( + this->fe_interface->n_current_interface_dofs()); + this->get_local_dof_values(fe_function, local_dof_values); + + get_jump_in_function_gradients_from_local_dof_values(local_dof_values, + gradients); + } + + + + template + template + void + Scalar::get_average_of_function_values_from_local_dof_values( + const InputVector &local_dof_values, + std::vector> &values) + const + { + Assert(values.size() == this->fe_interface->n_quadrature_points, + ExcNotImplemented()); + + const auto &interface_dof_indices = + this->fe_interface->get_interface_dof_indices(); + + for (const auto dof_index : this->fe_interface->dof_indices()) + for (const auto q_index : this->fe_interface->quadrature_point_indices()) + { + if (dof_index == 0) + values[q_index] = 0.; + + values[q_index] += + local_dof_values[dof_index] * average_of_values(dof_index, q_index); + } + } + + + + template + template + void + Scalar::get_average_of_function_values( + const InputVector &fe_function, + std::vector> &values) + const + { + std::vector local_dof_values( + this->fe_interface->n_current_interface_dofs()); + this->get_local_dof_values(fe_function, local_dof_values); + + get_average_of_function_values_from_local_dof_values(local_dof_values, + values); + } + + + + template + template + void + Scalar:: + get_average_of_function_gradients_from_local_dof_values( + const InputVector &local_dof_values, + std::vector> + &gradients) const + { + Assert(gradients.size() == this->fe_interface->n_quadrature_points, + ExcNotImplemented()); + + const auto &interface_dof_indices = + this->fe_interface->get_interface_dof_indices(); + + for (const auto dof_index : this->fe_interface->dof_indices()) + for (const auto q_index : this->fe_interface->quadrature_point_indices()) + { + if (dof_index == 0) + gradients[q_index] = 0.; + + gradients[q_index] += local_dof_values[dof_index] * + average_of_gradients(dof_index, q_index); + } + } + + + + template + template + void + Scalar::get_average_of_function_gradients( + const InputVector &fe_function, + std::vector> + &gradients) const + { + std::vector local_dof_values( + this->fe_interface->n_current_interface_dofs()); + this->get_local_dof_values(fe_function, local_dof_values); + + get_average_of_function_gradients_from_local_dof_values(local_dof_values, + gradients); + } + + + + template + template + void + Scalar::get_jump_in_function_hessians_from_local_dof_values( + const InputVector &local_dof_values, + std::vector> + &hessians) const + { + Assert(hessians.size() == this->fe_interface->n_quadrature_points, + ExcNotImplemented()); + + const auto &interface_dof_indices = + this->fe_interface->get_interface_dof_indices(); + + for (const auto dof_index : this->fe_interface->dof_indices()) + for (const auto q_index : this->fe_interface->quadrature_point_indices()) + { + if (dof_index == 0) + hessians[q_index] = 0.; + + hessians[q_index] += + local_dof_values[dof_index] * jump_in_hessians(dof_index, q_index); + } + } + + + + template + template + void + Scalar::get_jump_in_function_hessians( + const InputVector &fe_function, + std::vector> + &hessians) const + { + std::vector local_dof_values( + this->fe_interface->n_current_interface_dofs()); + this->get_local_dof_values(fe_function, local_dof_values); + + get_jump_in_function_hessians_from_local_dof_values(local_dof_values, + hessians); + } + + + + template + template + void + Scalar::get_average_of_function_hessians_from_local_dof_values( + const InputVector &local_dof_values, + std::vector> + &hessians) const + { + Assert(hessians.size() == this->fe_interface->n_quadrature_points, + ExcNotImplemented()); + + const auto &interface_dof_indices = + this->fe_interface->get_interface_dof_indices(); + + for (unsigned int dof_index = 0; dof_index < interface_dof_indices.size(); + ++dof_index) + for (const auto q_index : this->fe_interface->quadrature_point_indices()) + { + if (dof_index == 0) + hessians[q_index] = 0.; + + hessians[q_index] += local_dof_values[dof_index] * + average_of_hessians(dof_index, q_index); + } + } + + + + template + template + void + Scalar::get_average_of_function_hessians( + const InputVector &fe_function, + std::vector> + &hessians) const + { + std::vector local_dof_values( + this->fe_interface->n_current_interface_dofs()); + this->get_local_dof_values(fe_function, local_dof_values); + + get_average_of_function_hessians_from_local_dof_values(local_dof_values, + hessians); + } + + + + template + template + void + Scalar:: + get_jump_in_function_third_derivatives_from_local_dof_values( + const InputVector &local_dof_values, + std::vector< + solution_third_derivative_type> + &third_derivatives) const + { + Assert(third_derivatives.size() == this->fe_interface->n_quadrature_points, + ExcNotImplemented()); + + const auto &interface_dof_indices = + this->fe_interface->get_interface_dof_indices(); + + for (unsigned int dof_index = 0; dof_index < interface_dof_indices.size(); + ++dof_index) + for (const auto q_index : this->fe_interface->quadrature_point_indices()) + { + if (dof_index == 0) + third_derivatives[q_index] = 0.; + + third_derivatives[q_index] += + local_dof_values[dof_index] * + jump_in_third_derivatives(dof_index, q_index); + } + } + + + + template + template + void + Scalar::get_jump_in_function_third_derivatives( + const InputVector &fe_function, + std::vector< + solution_third_derivative_type> + &third_derivatives) const + { + std::vector local_dof_values( + this->fe_interface->n_current_interface_dofs()); + this->get_local_dof_values(fe_function, local_dof_values); + + get_jump_in_function_third_derivatives_from_local_dof_values( + local_dof_values, third_derivatives); + } + + + template Vector::Vector( const FEInterfaceValues &fe_interface, @@ -2423,6 +3494,375 @@ namespace FEInterfaceViews { return jump_in_third_derivatives(interface_dof_index, q_point); } + + + + template + template + void + Vector::get_function_values_from_local_dof_values( + const bool here_or_there, + const InputVector &local_dof_values, + std::vector> &values) + const + { + Assert(values.size() == this->fe_interface->n_quadrature_points, + ExcNotImplemented()); + + const auto &interface_dof_indices = + this->fe_interface->get_interface_dof_indices(); + + for (const auto dof_index : this->fe_interface->dof_indices()) + for (const auto q_index : this->fe_interface->quadrature_point_indices()) + { + if (dof_index == 0) + values[q_index] = 0.; + + values[q_index] += local_dof_values[dof_index] * + value(here_or_there, dof_index, q_index); + } + } + + + + template + template + void + Vector::get_function_values( + const bool here_or_there, + const InputVector &fe_function, + std::vector> &values) + const + { + std::vector local_dof_values( + this->fe_interface->n_current_interface_dofs()); + this->get_local_dof_values(fe_function, local_dof_values); + + get_function_values_from_local_dof_values(here_or_there, + local_dof_values, + values); + } + + + + template + template + void + Vector::get_jump_in_function_values_from_local_dof_values( + const InputVector &local_dof_values, + std::vector> &values) + const + { + Assert(values.size() == this->fe_interface->n_quadrature_points, + ExcNotImplemented()); + + const auto &interface_dof_indices = + this->fe_interface->get_interface_dof_indices(); + + for (const auto dof_index : this->fe_interface->dof_indices()) + for (const auto q_index : this->fe_interface->quadrature_point_indices()) + { + if (dof_index == 0) + values[q_index] = 0.; + + values[q_index] += + local_dof_values[dof_index] * jump_in_values(dof_index, q_index); + } + } + + + + template + template + void + Vector::get_jump_in_function_values( + const InputVector &fe_function, + std::vector> &values) + const + { + std::vector local_dof_values( + this->fe_interface->n_current_interface_dofs()); + this->get_local_dof_values(fe_function, local_dof_values); + + get_jump_in_function_values_from_local_dof_values(local_dof_values, values); + } + + + + template + template + void + Vector::get_jump_in_function_gradients_from_local_dof_values( + const InputVector &local_dof_values, + std::vector> + &gradients) const + { + Assert(gradients.size() == this->fe_interface->n_quadrature_points, + ExcNotImplemented()); + + const auto &interface_dof_indices = + this->fe_interface->get_interface_dof_indices(); + + for (const auto dof_index : this->fe_interface->dof_indices()) + for (const auto q_index : this->fe_interface->quadrature_point_indices()) + { + if (dof_index == 0) + gradients[q_index] = 0.; + + gradients[q_index] += + local_dof_values[dof_index] * jump_in_gradients(dof_index, q_index); + } + } + + + + template + template + void + Vector::get_jump_in_function_gradients( + const InputVector &fe_function, + std::vector> + &gradients) const + { + std::vector local_dof_values( + this->fe_interface->n_current_interface_dofs()); + this->get_local_dof_values(fe_function, local_dof_values); + + get_jump_in_function_gradients_from_local_dof_values(local_dof_values, + gradients); + } + + + + template + template + void + Vector::get_average_of_function_values_from_local_dof_values( + const InputVector &local_dof_values, + std::vector> &values) + const + { + Assert(values.size() == this->fe_interface->n_quadrature_points, + ExcNotImplemented()); + + const auto &interface_dof_indices = + this->fe_interface->get_interface_dof_indices(); + + for (const auto dof_index : this->fe_interface->dof_indices()) + for (const auto q_index : this->fe_interface->quadrature_point_indices()) + { + if (dof_index == 0) + values[q_index] = 0.; + + values[q_index] += + local_dof_values[dof_index] * average_of_values(dof_index, q_index); + } + } + + + + template + template + void + Vector::get_average_of_function_values( + const InputVector &fe_function, + std::vector> &values) + const + { + std::vector local_dof_values( + this->fe_interface->n_current_interface_dofs()); + this->get_local_dof_values(fe_function, local_dof_values); + + get_average_of_function_values_from_local_dof_values(local_dof_values, + values); + } + + + + template + template + void + Vector:: + get_average_of_function_gradients_from_local_dof_values( + const InputVector &local_dof_values, + std::vector> + &gradients) const + { + Assert(gradients.size() == this->fe_interface->n_quadrature_points, + ExcNotImplemented()); + + const auto &interface_dof_indices = + this->fe_interface->get_interface_dof_indices(); + + for (const auto dof_index : this->fe_interface->dof_indices()) + for (const auto q_index : this->fe_interface->quadrature_point_indices()) + { + if (dof_index == 0) + gradients[q_index] = 0.; + + gradients[q_index] += local_dof_values[dof_index] * + average_of_gradients(dof_index, q_index); + } + } + + + + template + template + void + Vector::get_average_of_function_gradients( + const InputVector &fe_function, + std::vector> + &gradients) const + { + std::vector local_dof_values( + this->fe_interface->n_current_interface_dofs()); + this->get_local_dof_values(fe_function, local_dof_values); + + get_average_of_function_gradients_from_local_dof_values(local_dof_values, + gradients); + } + + + + template + template + void + Vector::get_jump_in_function_hessians_from_local_dof_values( + const InputVector &local_dof_values, + std::vector> + &hessians) const + { + Assert(hessians.size() == this->fe_interface->n_quadrature_points, + ExcNotImplemented()); + + const auto &interface_dof_indices = + this->fe_interface->get_interface_dof_indices(); + + for (const auto dof_index : this->fe_interface->dof_indices()) + for (const auto q_index : this->fe_interface->quadrature_point_indices()) + { + if (dof_index == 0) + hessians[q_index] = 0.; + + hessians[q_index] += + local_dof_values[dof_index] * jump_in_hessians(dof_index, q_index); + } + } + + + + template + template + void + Vector::get_jump_in_function_hessians( + const InputVector &fe_function, + std::vector> + &hessians) const + { + std::vector local_dof_values( + this->fe_interface->n_current_interface_dofs()); + this->get_local_dof_values(fe_function, local_dof_values); + + get_jump_in_function_hessians_from_local_dof_values(local_dof_values, + hessians); + } + + + + template + template + void + Vector::get_average_of_function_hessians_from_local_dof_values( + const InputVector &local_dof_values, + std::vector> + &hessians) const + { + Assert(hessians.size() == this->fe_interface->n_quadrature_points, + ExcNotImplemented()); + + const auto &interface_dof_indices = + this->fe_interface->get_interface_dof_indices(); + + for (unsigned int dof_index = 0; dof_index < interface_dof_indices.size(); + ++dof_index) + for (const auto q_index : this->fe_interface->quadrature_point_indices()) + { + if (dof_index == 0) + hessians[q_index] = 0.; + + hessians[q_index] += local_dof_values[dof_index] * + average_of_hessians(dof_index, q_index); + } + } + + + + template + template + void + Vector::get_average_of_function_hessians( + const InputVector &fe_function, + std::vector> + &hessians) const + { + std::vector local_dof_values( + this->fe_interface->n_current_interface_dofs()); + this->get_local_dof_values(fe_function, local_dof_values); + + get_average_of_function_hessians_from_local_dof_values(local_dof_values, + hessians); + } + + + + template + template + void + Vector:: + get_jump_in_function_third_derivatives_from_local_dof_values( + const InputVector &local_dof_values, + std::vector< + solution_third_derivative_type> + &third_derivatives) const + { + Assert(third_derivatives.size() == this->fe_interface->n_quadrature_points, + ExcNotImplemented()); + + const auto &interface_dof_indices = + this->fe_interface->get_interface_dof_indices(); + + for (unsigned int dof_index = 0; dof_index < interface_dof_indices.size(); + ++dof_index) + for (const auto q_index : this->fe_interface->quadrature_point_indices()) + { + if (dof_index == 0) + third_derivatives[q_index] = 0.; + + third_derivatives[q_index] += + local_dof_values[dof_index] * + jump_in_third_derivatives(dof_index, q_index); + } + } + + + + template + template + void + Vector::get_jump_in_function_third_derivatives( + const InputVector &fe_function, + std::vector< + solution_third_derivative_type> + &third_derivatives) const + { + std::vector local_dof_values( + this->fe_interface->n_current_interface_dofs()); + this->get_local_dof_values(fe_function, local_dof_values); + + get_jump_in_function_third_derivatives_from_local_dof_values( + local_dof_values, third_derivatives); + } } // namespace FEInterfaceViews #endif // DOXYGEN diff --git a/tests/feinterface/fe_interface_view_get_function_from_local_dof_values_01.cc b/tests/feinterface/fe_interface_view_get_function_from_local_dof_values_01.cc new file mode 100644 index 0000000000..21316bfd54 --- /dev/null +++ b/tests/feinterface/fe_interface_view_get_function_from_local_dof_values_01.cc @@ -0,0 +1,467 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2017 - 2021 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.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- + + +// Check that FEInterfaceViews::get_function_*_from_local_dof_values +// works with different number types. + +#include + +#include +#include + +#include +#include +#include +#include +#include + +#include +#include + +#include + +#include +#include + +#include "../tests.h" + +template +void +test_feiv_view(const Vector & solution, + const FEInterfaceValues & fe_iv, + const unsigned int & n_q_points, + const ExtractorType & extractor, + const std::vector &local_dof_values); + +// Scalar view +template +void +test_feiv_view(const Vector & solution, + const FEInterfaceValues & fe_iv, + const unsigned int & n_q_points, + const FEValuesExtractors::Scalar &extractor, + const std::vector & local_dof_values) +{ + using View = typename std::remove_reference< + typename std::remove_const::type>::type; + const View &fe_iv_view = fe_iv[extractor]; + + // Typedefs + using value_type = + typename ProductType::type; + using gradient_type = + typename ProductType::type; + using hessian_type = + typename ProductType::type; + using laplacian_type = + typename ProductType::type; + using third_derivative_type = + typename ProductType::type; + + // Values + std::vector> + qp_jumps_local(n_q_points), qp_averages_local(n_q_points), + qp_interface_values_local(n_q_points); + std::vector qp_jumps_global(n_q_points), + qp_averages_global(n_q_points), qp_interface_values_global(n_q_points); + + fe_iv_view.get_jump_in_function_values_from_local_dof_values(local_dof_values, + qp_jumps_local); + fe_iv_view.get_jump_in_function_values(solution, qp_jumps_global); + + fe_iv_view.get_average_of_function_values_from_local_dof_values( + local_dof_values, qp_averages_local); + fe_iv_view.get_average_of_function_values(solution, qp_averages_global); + + fe_iv_view.get_function_values_from_local_dof_values( + true, local_dof_values, qp_interface_values_local); + + fe_iv_view.get_function_values(true, solution, qp_interface_values_global); + + // Gradients + std::vector> + qp_jump_grads_local(n_q_points), qp_average_grads_local(n_q_points); + std::vector qp_jump_grads_global(n_q_points), + qp_average_grads_global(n_q_points); + + fe_iv_view.get_jump_in_function_gradients_from_local_dof_values( + local_dof_values, qp_jump_grads_local); + fe_iv_view.get_jump_in_function_gradients(solution, qp_jump_grads_global); + + fe_iv_view.get_average_of_function_gradients_from_local_dof_values( + local_dof_values, qp_average_grads_local); + fe_iv_view.get_average_of_function_gradients(solution, + qp_average_grads_global); + + // Hessians + std::vector> + qp_jump_hess_local(n_q_points), qp_average_hess_local(n_q_points); + std::vector qp_jump_hess_global(n_q_points), + qp_average_hess_global(n_q_points); + + fe_iv_view.get_jump_in_function_hessians_from_local_dof_values( + local_dof_values, qp_jump_hess_local); + fe_iv_view.get_jump_in_function_hessians(solution, qp_jump_hess_global); + fe_iv_view.get_average_of_function_hessians_from_local_dof_values( + local_dof_values, qp_average_hess_local); + fe_iv_view.get_average_of_function_hessians(solution, qp_average_hess_global); + + // Third derivatives + std::vector< + typename View::template solution_third_derivative_type> + qp_jump_third_deriv_local(n_q_points); + std::vector qp_jump_third_deriv_global(n_q_points); + + fe_iv_view.get_jump_in_function_third_derivatives_from_local_dof_values( + local_dof_values, qp_jump_third_deriv_local); + fe_iv_view.get_jump_in_function_third_derivatives(solution, + qp_jump_third_deriv_global); + + + // Output + for (unsigned int q = 0; q < n_q_points; ++q) + { + if (value_type(qp_interface_values_local[q]) != + value_type(qp_interface_values_global[q])) + deallog << "Interface value NOT OK: Value @ " << q << std::endl; + + if (value_type(qp_jumps_local[q]) != value_type(qp_jumps_global[q])) + deallog << "Jump NOT OK: Value @ " << q << std::endl; + + if (value_type(qp_averages_local[q]) != value_type(qp_averages_global[q])) + deallog << "Average NOT OK: Value @ " << q << std::endl; + + if (gradient_type(qp_jump_grads_local[q]) != + gradient_type(qp_jump_grads_global[q])) + deallog << "Jump NOT OK: Grad @ " << q << std::endl; + + if (gradient_type(qp_average_grads_local[q]) != + gradient_type(qp_average_grads_global[q])) + deallog << "Average NOT OK: Grad @ " << q << std::endl; + + if (hessian_type(qp_jump_hess_local[q]) != + hessian_type(qp_jump_hess_global[q])) + deallog << "Jump NOT OK: Hess @ " << q << std::endl; + + if (hessian_type(qp_average_hess_local[q]) != + hessian_type(qp_average_hess_global[q])) + deallog << "Average NOT OK: Hess @ " << q << std::endl; + + if (third_derivative_type(qp_jump_third_deriv_local[q]) != + third_derivative_type(qp_jump_third_deriv_global[q])) + deallog << "Jump NOT OK: 3rd der @ " << q << std::endl; + } +} + +// Vector view +template +void +test_feiv_view(const Vector & solution, + const FEInterfaceValues & fe_iv, + const unsigned int & n_q_points, + const FEValuesExtractors::Vector &extractor, + const std::vector & local_dof_values) +{ + using View = typename std::remove_reference< + typename std::remove_const::type>::type; + const View &fe_iv_view = fe_iv[extractor]; + + // Typedefs + using value_type = + typename ProductType::type; + using gradient_type = + typename ProductType::type; + using hessian_type = + typename ProductType::type; + using third_derivative_type = + typename ProductType::type; + + // Values + std::vector> + qp_jumps_local(n_q_points), qp_averages_local(n_q_points), + qp_interface_values_local(n_q_points); + std::vector qp_jumps_global(n_q_points), + qp_averages_global(n_q_points), qp_interface_values_global(n_q_points); + + fe_iv_view.get_jump_in_function_values_from_local_dof_values(local_dof_values, + qp_jumps_local); + fe_iv_view.get_jump_in_function_values(solution, qp_jumps_global); + + fe_iv_view.get_average_of_function_values_from_local_dof_values( + local_dof_values, qp_averages_local); + fe_iv_view.get_average_of_function_values(solution, qp_averages_global); + + fe_iv_view.get_function_values_from_local_dof_values( + true, local_dof_values, qp_interface_values_local); + + fe_iv_view.get_function_values(true, solution, qp_interface_values_global); + + // Gradients + std::vector> + qp_jump_grads_local(n_q_points), qp_average_grads_local(n_q_points); + std::vector qp_jump_grads_global(n_q_points), + qp_average_grads_global(n_q_points); + + fe_iv_view.get_jump_in_function_gradients_from_local_dof_values( + local_dof_values, qp_jump_grads_local); + fe_iv_view.get_jump_in_function_gradients(solution, qp_jump_grads_global); + + fe_iv_view.get_average_of_function_gradients_from_local_dof_values( + local_dof_values, qp_average_grads_local); + fe_iv_view.get_average_of_function_gradients(solution, + qp_average_grads_global); + + // Hessians + std::vector> + qp_jump_hess_local(n_q_points), qp_average_hess_local(n_q_points); + std::vector qp_jump_hess_global(n_q_points), + qp_average_hess_global(n_q_points); + + fe_iv_view.get_jump_in_function_hessians_from_local_dof_values( + local_dof_values, qp_jump_hess_local); + fe_iv_view.get_jump_in_function_hessians(solution, qp_jump_hess_global); + fe_iv_view.get_average_of_function_hessians_from_local_dof_values( + local_dof_values, qp_average_hess_local); + fe_iv_view.get_average_of_function_hessians(solution, qp_average_hess_global); + + // Third derivatives + std::vector< + typename View::template solution_third_derivative_type> + qp_jump_third_deriv_local(n_q_points); + std::vector qp_jump_third_deriv_global(n_q_points); + + fe_iv_view.get_jump_in_function_third_derivatives_from_local_dof_values( + local_dof_values, qp_jump_third_deriv_local); + fe_iv_view.get_jump_in_function_third_derivatives(solution, + qp_jump_third_deriv_global); + + + // Output + for (unsigned int q = 0; q < n_q_points; ++q) + { + if (value_type(qp_interface_values_local[q]) != + value_type(qp_interface_values_global[q])) + deallog << "Interface value NOT OK: Value @ " << q << std::endl; + + if (value_type(qp_jumps_local[q]) != value_type(qp_jumps_global[q])) + deallog << "Jump NOT OK: Value @ " << q << std::endl; + + if (value_type(qp_averages_local[q]) != value_type(qp_averages_global[q])) + deallog << "Average NOT OK: Value @ " << q << std::endl; + + if (value_type(qp_averages_local[q]) != value_type(qp_averages_global[q])) + deallog << "Average NOT OK: Value @ " << q << std::endl; + + if (gradient_type(qp_jump_grads_local[q]) != + gradient_type(qp_jump_grads_global[q])) + deallog << "Jump NOT OK: Grad @ " << q << std::endl; + + if (gradient_type(qp_average_grads_local[q]) != + gradient_type(qp_average_grads_global[q])) + deallog << "Average NOT OK: Grad @ " << q << std::endl; + + if (hessian_type(qp_jump_hess_local[q]) != + hessian_type(qp_jump_hess_global[q])) + deallog << "Jump NOT OK: Hess @ " << q << std::endl; + + if (hessian_type(qp_average_hess_local[q]) != + hessian_type(qp_average_hess_global[q])) + deallog << "Average NOT OK: Hess @ " << q << std::endl; + + if (third_derivative_type(qp_jump_third_deriv_local[q]) != + third_derivative_type(qp_jump_third_deriv_global[q])) + deallog << "Jump NOT OK: 3rd der @ " << q << std::endl; + } +} + +template +void +test_feiv_extractor(const FEType &fe, const ExtractorType &extractor) +{ + Triangulation tria; + DoFHandler dof_handler(tria); + Vector solution; + + GridGenerator::hyper_cube(tria, -1, 1); + tria.refine_global(1); + dof_handler.distribute_dofs(fe); + DoFRenumbering::random(dof_handler); + solution.reinit(dof_handler.n_dofs()); + + // Populate with non-trivial values + for (unsigned int i = 0; i < dof_handler.n_dofs(); ++i) + { + solution(i) = i + 1; + } + + const QGauss face_quadrature(2); + FEInterfaceValues fe_iv(fe, + face_quadrature, + update_values | update_gradients | + update_quadrature_points | + update_3rd_derivatives); + + typename DoFHandler::active_cell_iterator cell = + dof_handler.begin_active(); + for (const auto f : cell->face_indices()) + if (!cell->face(f)->at_boundary()) + { + { + fe_iv.reinit(cell, + f, + numbers::invalid_unsigned_int, + cell->neighbor(f), + cell->neighbor_of_neighbor(f), + numbers::invalid_unsigned_int); + const std::vector &interface_dof_indices = + fe_iv.get_interface_dof_indices(); + + // Convert the DoF values so that they are potentially of + // a different number type + std::vector interface_dof_values_other( + fe_iv.n_current_interface_dofs()); + for (unsigned int i = 0; i < fe_iv.n_current_interface_dofs(); ++i) + interface_dof_values_other[i] = solution(interface_dof_indices[i]); + + test_feiv_view(solution, + fe_iv, + face_quadrature.size(), + extractor, + interface_dof_values_other); + } + } + + deallog << "OK" << std::endl; +} + + +template +void +test_jump_function() +{ + FE_DGP fe(0); + FEValuesExtractors::Scalar extractor(0); + + Triangulation tria; + DoFHandler dof_handler(tria); + Vector solution; + + GridGenerator::hyper_cube(tria, -1, 1); + tria.refine_global(1); + dof_handler.distribute_dofs(fe); + + solution.reinit(dof_handler.n_dofs()); + + std::vector local_dof_indices(1); + // Populate with non-trivial values + for (const auto &cell : dof_handler.active_cell_iterators()) + { + cell->get_dof_indices(local_dof_indices); + solution(local_dof_indices[0]) = + static_cast(cell->active_cell_index()); + } + + const QGauss face_quadrature(2); + FEInterfaceValues fe_iv(fe, + face_quadrature, + update_values | update_gradients | + update_quadrature_points | + update_3rd_derivatives); + + std::vector qp_jumps_global(face_quadrature.size()); + + typename DoFHandler::active_cell_iterator cell = + dof_handler.begin_active(); + for (; cell != dof_handler.end(); ++cell) + for (const auto f : cell->face_indices()) + if (!cell->face(f)->at_boundary()) + { + { + fe_iv.reinit(cell, + f, + numbers::invalid_unsigned_int, + cell->neighbor(f), + cell->neighbor_of_neighbor(f), + numbers::invalid_unsigned_int); + + fe_iv[extractor].get_jump_in_function_values(solution, + qp_jumps_global); + + double exact = cell->active_cell_index() * 1.0 - + cell->neighbor(f)->active_cell_index() * 1.0; + + for (unsigned int q = 0; q < face_quadrature.size(); ++q) + Assert(std::fabs(qp_jumps_global[q] - exact) < 1e-15, + ExcNotImplemented()); + } + } + + deallog << "OK" << std::endl; +} + +template +void +test() +{ + const unsigned int degree = 3; // Need third derivatives + + deallog.push("FEInterfaceViews Scalar"); + { + FE_Q fe(degree); + FEValuesExtractors::Scalar extractor(0); + test_feiv_extractor(fe, extractor); + } + deallog.pop(); + + deallog.push("FEInterfaceViews Vector"); + { + FESystem fe(FE_Q(degree), dim); + FEValuesExtractors::Vector extractor(0); + test_feiv_extractor(fe, extractor); + } + deallog.pop(); +} + +int +main(int argc, char **argv) +{ + initlog(); + + Utilities::MPI::MPI_InitFinalize mpi_initialization( + argc, argv, testing_max_num_threads()); + + deallog.push("Float"); + { + test(); + } + deallog.pop(); + + deallog.push("Double"); + { + test(); + } + deallog.pop(); + + deallog.push("Test jump functions"); + { + test_jump_function(); + } + deallog.pop(); + + deallog << "OK" << std::endl; +} diff --git a/tests/feinterface/fe_interface_view_get_function_from_local_dof_values_01.output b/tests/feinterface/fe_interface_view_get_function_from_local_dof_values_01.output new file mode 100644 index 0000000000..67e69f6ac8 --- /dev/null +++ b/tests/feinterface/fe_interface_view_get_function_from_local_dof_values_01.output @@ -0,0 +1,7 @@ + +DEAL:Float:FEInterfaceViews Scalar::OK +DEAL:Float:FEInterfaceViews Vector::OK +DEAL:Double:FEInterfaceViews Scalar::OK +DEAL:Double:FEInterfaceViews Vector::OK +DEAL:Test jump functions::OK +DEAL::OK diff --git a/tests/feinterface/fe_interface_view_get_function_from_local_dof_values_02.cc b/tests/feinterface/fe_interface_view_get_function_from_local_dof_values_02.cc new file mode 100644 index 0000000000..d54d01389b --- /dev/null +++ b/tests/feinterface/fe_interface_view_get_function_from_local_dof_values_02.cc @@ -0,0 +1,321 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2017 - 2021 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.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- + + + +// Check that FEInterfaceViews::get_function_*_from_local_dof_values +// works with different number types. +// As opposed to fe_interface_view_get_function_from_local_dof_values_01.cc, +// this test does not check the output types, but rather that all necessary +// instantiations are in place and that no assertions are triggered. +// This is because FEInterfaceViews::get_function_* cannot work with +// Sacado types. + +#include + +#include + +#include +#include + +#include +#include +#include +#include + +#include +#include + +#include + +#include +#include + +#include "../tests.h" + + +template +void +test_feiv_view(const Vector & solution, + const FEValues & fe_values, + const unsigned int & n_q_points, + const ExtractorType & extractor, + const std::vector &local_dof_values); + +// Scalar view +template +void +test_feiv_view(const Vector & solution, + const FEInterfaceValues & fe_iv, + const unsigned int & n_q_points, + const FEValuesExtractors::Scalar &extractor, + const std::vector & local_dof_values) +{ + using View = typename std::remove_reference< + typename std::remove_const::type>::type; + const View &fe_iv_view = fe_iv[extractor]; + + // Typedefs + using value_type = + typename ProductType::type; + using gradient_type = + typename ProductType::type; + using hessian_type = + typename ProductType::type; + using third_derivative_type = + typename ProductType::type; + + // Values + std::vector> + qp_jumps_local(n_q_points), qp_averages_local(n_q_points), + qp_interface_values_local(n_q_points); + std::vector qp_jumps_global(n_q_points), + qp_averages_global(n_q_points), qp_interface_values_global(n_q_points); + + fe_iv_view.get_jump_in_function_values_from_local_dof_values(local_dof_values, + qp_jumps_local); + + fe_iv_view.get_average_of_function_values_from_local_dof_values( + local_dof_values, qp_averages_local); + + fe_iv_view.get_function_values_from_local_dof_values( + true, local_dof_values, qp_interface_values_local); + + // Gradients + std::vector> + qp_jump_grads_local(n_q_points), qp_average_grads_local(n_q_points); + std::vector qp_jump_grads_global(n_q_points), + qp_average_grads_global(n_q_points); + + fe_iv_view.get_jump_in_function_gradients_from_local_dof_values( + local_dof_values, qp_jump_grads_local); + + fe_iv_view.get_average_of_function_gradients_from_local_dof_values( + local_dof_values, qp_average_grads_local); + + // Hessians + std::vector> + qp_jump_hess_local(n_q_points), qp_average_hess_local(n_q_points); + std::vector qp_jump_hess_global(n_q_points), + qp_average_hess_global(n_q_points); + + fe_iv_view.get_jump_in_function_hessians_from_local_dof_values( + local_dof_values, qp_jump_hess_local); + + fe_iv_view.get_average_of_function_hessians_from_local_dof_values( + local_dof_values, qp_average_hess_local); + + // Third derivatives + std::vector< + typename View::template solution_third_derivative_type> + qp_jump_third_deriv_local(n_q_points); + std::vector qp_jump_third_deriv_global(n_q_points); + + fe_iv_view.get_jump_in_function_third_derivatives_from_local_dof_values( + local_dof_values, qp_jump_third_deriv_local); +} + +// Vector view +template +void +test_feiv_view(const Vector & solution, + const FEInterfaceValues & fe_iv, + const unsigned int & n_q_points, + const FEValuesExtractors::Vector &extractor, + const std::vector & local_dof_values) +{ + using View = typename std::remove_reference< + typename std::remove_const::type>::type; + const View &fe_iv_view = fe_iv[extractor]; + + // Typedefs + using value_type = + typename ProductType::type; + using gradient_type = + typename ProductType::type; + using hessian_type = + typename ProductType::type; + using third_derivative_type = + typename ProductType::type; + + // Values + std::vector> + qp_jumps_local(n_q_points), qp_averages_local(n_q_points), + qp_interface_values_local(n_q_points); + std::vector qp_jumps_global(n_q_points), + qp_averages_global(n_q_points), qp_interface_values_global(n_q_points); + + fe_iv_view.get_jump_in_function_values_from_local_dof_values(local_dof_values, + qp_jumps_local); + + fe_iv_view.get_average_of_function_values_from_local_dof_values( + local_dof_values, qp_averages_local); + + fe_iv_view.get_function_values_from_local_dof_values( + true, local_dof_values, qp_interface_values_local); + + // Gradients + std::vector> + qp_jump_grads_local(n_q_points), qp_average_grads_local(n_q_points); + std::vector qp_jump_grads_global(n_q_points), + qp_average_grads_global(n_q_points); + + fe_iv_view.get_jump_in_function_gradients_from_local_dof_values( + local_dof_values, qp_jump_grads_local); + + fe_iv_view.get_average_of_function_gradients_from_local_dof_values( + local_dof_values, qp_average_grads_local); + + // Hessians + std::vector> + qp_jump_hess_local(n_q_points), qp_average_hess_local(n_q_points); + std::vector qp_jump_hess_global(n_q_points), + qp_average_hess_global(n_q_points); + + fe_iv_view.get_jump_in_function_hessians_from_local_dof_values( + local_dof_values, qp_jump_hess_local); + + fe_iv_view.get_average_of_function_hessians_from_local_dof_values( + local_dof_values, qp_average_hess_local); + + // Third derivatives + std::vector< + typename View::template solution_third_derivative_type> + qp_jump_third_deriv_local(n_q_points); + std::vector qp_jump_third_deriv_global(n_q_points); + + fe_iv_view.get_jump_in_function_third_derivatives_from_local_dof_values( + local_dof_values, qp_jump_third_deriv_local); +} + +template +void +test_feiv_extractor(const FEType &fe, const ExtractorType &extractor) +{ + Triangulation tria; + DoFHandler dof_handler(tria); + Vector solution; + + GridGenerator::hyper_cube(tria, -1, 1); + tria.refine_global(1); + dof_handler.distribute_dofs(fe); + DoFRenumbering::random(dof_handler); + solution.reinit(dof_handler.n_dofs()); + + // Populate with non-trivial values + for (unsigned int i = 0; i < dof_handler.n_dofs(); ++i) + { + solution(i) = i + 1; + } + + const QGauss face_quadrature(2); + FEInterfaceValues fe_iv(fe, + face_quadrature, + update_values | update_gradients | + update_quadrature_points | + update_3rd_derivatives); + + typename DoFHandler::active_cell_iterator cell = + dof_handler.begin_active(); + for (const auto f : cell->face_indices()) + if (!cell->face(f)->at_boundary()) + { + { + fe_iv.reinit(cell, + f, + numbers::invalid_unsigned_int, + cell->neighbor(f), + cell->neighbor_of_neighbor(f), + numbers::invalid_unsigned_int); + const std::vector &interface_dof_indices = + fe_iv.get_interface_dof_indices(); + + // Convert the DoF values so that they are potentially of + // a different number type + std::vector interface_dof_values_other( + fe_iv.n_current_interface_dofs()); + for (unsigned int i = 0; i < fe_iv.n_current_interface_dofs(); ++i) + interface_dof_values_other[i] = solution(interface_dof_indices[i]); + + test_feiv_view(solution, + fe_iv, + face_quadrature.size(), + extractor, + interface_dof_values_other); + } + } + + deallog << "OK" << std::endl; +} + +template +void +test() +{ + const unsigned int degree = 3; // Need third derivatives + + deallog.push("FEInterfaceViews Scalar"); + { + FE_Q fe(degree); + FEValuesExtractors::Scalar extractor(0); + test_feiv_extractor(fe, extractor); + } + deallog.pop(); + + deallog.push("FEInterfaceViews Vector"); + { + FESystem fe(FE_Q(degree), dim); + FEValuesExtractors::Vector extractor(0); + test_feiv_extractor(fe, extractor); + } + deallog.pop(); +} + +int +main(int argc, char **argv) +{ + initlog(); + + Utilities::MPI::MPI_InitFinalize mpi_initialization( + argc, argv, testing_max_num_threads()); + + deallog.push("Sacado::Fad::DFad"); + { + test>(); + } + deallog.pop(); + + deallog.push("Sacado::Fad::DFad>"); + { + test>>(); + } + deallog.pop(); + + deallog.push("Sacado::Fad::DFad"); + { + test>(); + } + deallog.pop(); + + deallog.push("Sacado::Fad::DFad>"); + { + test>>(); + } + deallog.pop(); + + deallog << "OK" << std::endl; +} diff --git a/tests/feinterface/fe_interface_view_get_function_from_local_dof_values_02.with_trilinos=true.with_trilinos_wiith_sacado=true.output b/tests/feinterface/fe_interface_view_get_function_from_local_dof_values_02.with_trilinos=true.with_trilinos_wiith_sacado=true.output new file mode 100644 index 0000000000..1fb2479956 --- /dev/null +++ b/tests/feinterface/fe_interface_view_get_function_from_local_dof_values_02.with_trilinos=true.with_trilinos_wiith_sacado=true.output @@ -0,0 +1,10 @@ + +DEAL:Sacado::Fad::DFad:FEInterfaceViews Scalar::OK +DEAL:Sacado::Fad::DFad:FEInterfaceViews Vector::OK +DEAL:Sacado::Fad::DFad>:FEInterfaceViews Scalar::OK +DEAL:Sacado::Fad::DFad>:FEInterfaceViews Vector::OK +DEAL:Sacado::Fad::DFad:FEInterfaceViews Scalar::OK +DEAL:Sacado::Fad::DFad:FEInterfaceViews Vector::OK +DEAL:Sacado::Fad::DFad>:FEInterfaceViews Scalar::OK +DEAL:Sacado::Fad::DFad>:FEInterfaceViews Vector::OK +DEAL::OK