From a81099a807b7dff82d29dbe43d7d492fd4e500c3 Mon Sep 17 00:00:00 2001 From: Jean-Paul Pelteret Date: Sat, 14 Aug 2021 18:38:31 +0200 Subject: [PATCH] Extend ScratchData to return jumps and averages of FE functions on an interface --- include/deal.II/fe/fe_interface_values.h | 52 +++ include/deal.II/meshworker/scratch_data.h | 536 +++++++++++++++++++++- source/meshworker/scratch_data.cc | 12 + 3 files changed, 592 insertions(+), 8 deletions(-) diff --git a/include/deal.II/fe/fe_interface_values.h b/include/deal.II/fe/fe_interface_values.h index d5b228d2d0..5447a074b3 100644 --- a/include/deal.II/fe/fe_interface_values.h +++ b/include/deal.II/fe/fe_interface_values.h @@ -1175,6 +1175,58 @@ namespace FEInterfaceViews } // namespace FEInterfaceViews +namespace internal +{ + namespace FEInterfaceViews + { + /** + * A class whose specialization is used to define what FEInterfaceViews + * object corresponds to the given FEValuesExtractors object. + */ + template + struct ViewType + {}; + + /** + * A class whose specialization is used to define what FEInterfaceViews + * object corresponds to the given FEValuesExtractors object. + * + * When using FEValuesExtractors::Scalar, the corresponding view is an + * FEInterfaceViews::Scalar. + */ + template + struct ViewType + { + using type = typename dealii::FEInterfaceViews::Scalar; + }; + + /** + * A class whose specialization is used to define what FEInterfaceViews + * object corresponds to the given FEValuesExtractors object. + * + * When using FEValuesExtractors::Vector, the corresponding view is an + * FEInterfaceViews::Vector. + */ + template + struct ViewType + { + using type = typename dealii::FEInterfaceViews::Vector; + }; + } // namespace FEInterfaceViews +} // namespace internal + +namespace FEInterfaceViews +{ + /** + * A templated alias that associates to a given Extractor class + * the corresponding view in FEInterfaceViews. + */ + template + using View = typename dealii::internal::FEInterfaceViews:: + ViewType::type; +} // namespace FEInterfaceViews + + /** * FEInterfaceValues is a data structure to access and assemble finite element diff --git a/include/deal.II/meshworker/scratch_data.h b/include/deal.II/meshworker/scratch_data.h index 318a64cf19..dc7bbd433d 100644 --- a/include/deal.II/meshworker/scratch_data.h +++ b/include/deal.II/meshworker/scratch_data.h @@ -372,7 +372,7 @@ namespace MeshWorker * get_current_fe_values() will return the FEValuesBase associated to the * current cell, while get_neighbor_fe_values() will be associated with the * neighbor cell. The method get_local_dof_indices() will return the - * same result of FEInterfaceValues::get_interface_dof_to_dof_indices(), + * same result of FEInterfaceValues::get_interface_dof_indices(), * while the get_neighbor_dof_indices() will return the local dof indices * of the neighbor cell. */ @@ -398,6 +398,12 @@ namespace MeshWorker const FEValuesBase & get_current_fe_values() const; + /** + * Get the currently initialized FEInterfaceValues. + */ + const FEInterfaceValues & + get_current_interface_fe_values() const; + /** * Return the quadrature points of the internal FEValues object. */ @@ -520,6 +526,12 @@ namespace MeshWorker const GeneralDataStorage & get_general_data_storage() const; + /** + * Return a reference to the used mapping. + */ + const Mapping & + get_mapping() const; + /** * @name Evaluation of finite element fields and their derivatives on the current cell */ @@ -716,7 +728,7 @@ namespace MeshWorker /** * For the solution vector identified by @p global_vector_name, compute - * the hessians of the function at the quadrature points, and return a + * the Hessians of the function at the quadrature points, and return a * vector with the correct type deduced by the Extractor you passed as the * @p variable argument. * @@ -772,7 +784,7 @@ namespace MeshWorker /** * For the solution vector identified by @p global_vector_name, compute - * the third_derivatives of the function at the quadrature points, and + * the third derivatives of the function at the quadrature points, and * return a vector with the correct type deduced by the Extractor you passed * as the @p variable argument. * @@ -801,10 +813,228 @@ namespace MeshWorker /** @} */ // CurrentCellEvaluation /** - * Return a reference to the used mapping. + * @name Evaluation of jumps in finite element fields and their derivatives on the current interface */ - const Mapping & - get_mapping() const; + /** @{ */ // CurrentInterfaceJumpEvaluation + + /** + * For the solution vector identified by @p global_vector_name, compute + * the jumps in the values of the function at the quadrature points, and + * return a vector with the correct type deduced by the interface Extractor + * you passed as the + * @p variable argument. + * + * Before you can call this method, you need to call the + * extract_local_dof_values() method at least once, passing the same + * @p global_vector_name string, and the same type for the variable @p dummy. + * + * If you have not previously called the extract_local_dof_values() method, + * this function will throw an exception. + * + * For this function to work properly, the underlying FEInterfaceValues + * object for which you called one of the reinit() functions, must have + * computed the information you are requesting. To do so, the update_values + * flag must be an element of the list of UpdateFlags that you passed to the + * constructor of this object. See "The interplay of UpdateFlags, Mapping, + * and FiniteElement" in the documentation of the FEValues class for more + * information. + */ + template + const std::vector< + typename FEInterfaceViews::View:: + template solution_value_type> & + get_jumps_in_values(const std::string &global_vector_name, + const Extractor & variable, + const Number dummy = Number(0)); + + /** + * For the solution vector identified by @p global_vector_name, compute + * the jumps in the gradients of the function at the quadrature points, and + * return a vector with the correct type deduced by the interface Extractor + * you passed as the + * @p variable argument. + * + * Before you can call this method, you need to call the + * extract_local_dof_values() method at least once, passing the same + * @p global_vector_name string, and the same type for the variable @p dummy. + * + * If you have not previously called the extract_local_dof_values() method, + * this function will throw an exception. + * + * For this function to work properly, the underlying FEInterfaceValues + * object for which you called one of the reinit() functions, must have + * computed the information you are requesting. To do so, the + * update_gradients flag must be an element of the list of UpdateFlags that + * you passed to the constructor of this object. See "The interplay of + * UpdateFlags, Mapping, and FiniteElement" in the documentation of the + * FEValues class for more information. + */ + template + const std::vector< + typename FEInterfaceViews::View:: + template solution_gradient_type> & + get_jumps_in_gradients(const std::string &global_vector_name, + const Extractor & variable, + const Number dummy = Number(0)); + + /** + * For the solution vector identified by @p global_vector_name, compute + * the jumps in the Hessians of the function at the quadrature points, and + * return a vector with the correct type deduced by the interface Extractor + * you passed as the + * @p variable argument. + * + * Before you can call this method, you need to call the + * extract_local_dof_values() method at least once, passing the same + * @p global_vector_name string, and the same type for the variable @p dummy. + * + * If you have not previously called the extract_local_dof_values() method, + * this function will throw an exception. + * + * For this function to work properly, the underlying FEInterfaceValues + * object for which you called one of the reinit() functions, must have + * computed the information you are requesting. To do so, the + * update_hessians flag must be an element of the list of UpdateFlags that + * you passed to the constructor of this object. See "The interplay of + * UpdateFlags, Mapping, and FiniteElement" in the documentation of the + * FEValues class for more information. + */ + template + const std::vector< + typename FEInterfaceViews::View:: + template solution_hessian_type> & + get_jumps_in_hessians(const std::string &global_vector_name, + const Extractor & variable, + const Number dummy = Number(0)); + + /** + * For the solution vector identified by @p global_vector_name, compute + * the jumps in the third derivatives of the function at the quadrature + * points, and return a vector with the correct type deduced by the + * interface Extractor you passed + * as the @p variable argument. + * + * Before you can call this method, you need to call the + * extract_local_dof_values() method at least once, passing the same + * @p global_vector_name string, and the same type for the variable @p dummy. + * + * If you have not previously called the extract_local_dof_values() method, + * this function will throw an exception. + * + * For this function to work properly, the underlying FEInterfaceValues + * object for which you called one of the reinit() functions, must have + * computed the information you are requesting. To do so, the + * update_3rd_derivatives flag must be an element of the list of UpdateFlags + * that you passed to the constructor of this object. See "The interplay of + * UpdateFlags, Mapping, and FiniteElement" in the documentation of the + * FEValues for more information. + */ + template + const std::vector< + typename FEInterfaceViews::View:: + template solution_third_derivative_type> & + get_jumps_in_third_derivatives(const std::string &global_vector_name, + const Extractor & variable, + const Number dummy = Number(0)); + + /** @} */ // CurrentInterfaceJumpEvaluation + + /** + * @name Evaluation of averages of finite element fields and their derivatives on the current interface + */ + /** @{ */ // CurrentInterfaceAverageEvaluation + + /** + * For the solution vector identified by @p global_vector_name, compute + * the averages of the values of the function at the quadrature points, and + * return a vector with the correct type deduced by the interface Extractor + * you passed as the + * @p variable argument. + * + * Before you can call this method, you need to call the + * extract_local_dof_values() method at least once, passing the same + * @p global_vector_name string, and the same type for the variable @p dummy. + * + * If you have not previously called the extract_local_dof_values() method, + * this function will throw an exception. + * + * For this function to work properly, the underlying FEInterfaceValues + * object for which you called one of the reinit() functions, must have + * computed the information you are requesting. To do so, the update_values + * flag must be an element of the list of UpdateFlags that you passed to the + * constructor of this object. See "The interplay of UpdateFlags, Mapping, + * and FiniteElement" in the documentation of the FEValues class for more + * information. + */ + template + const std::vector< + typename FEInterfaceViews::View:: + template solution_value_type> & + get_averages_of_values(const std::string &global_vector_name, + const Extractor & variable, + const Number dummy = Number(0)); + + /** + * For the solution vector identified by @p global_vector_name, compute + * the averages of the gradients of the function at the quadrature points, + * and return a vector with the correct type deduced by the interface + * Extractor you passed as the + * @p variable argument. + * + * Before you can call this method, you need to call the + * extract_local_dof_values() method at least once, passing the same + * @p global_vector_name string, and the same type for the variable @p dummy. + * + * If you have not previously called the extract_local_dof_values() method, + * this function will throw an exception. + * + * For this function to work properly, the underlying FEInterfaceValues + * object for which you called one of the reinit() functions, must have + * computed the information you are requesting. To do so, the + * update_gradients flag must be an element of the list of UpdateFlags that + * you passed to the constructor of this object. See "The interplay of + * UpdateFlags, Mapping, and FiniteElement" in the documentation of the + * FEValues class for more information. + */ + template + const std::vector< + typename FEInterfaceViews::View:: + template solution_gradient_type> & + get_averages_of_gradients(const std::string &global_vector_name, + const Extractor & variable, + const Number dummy = Number(0)); + + /** + * For the solution vector identified by @p global_vector_name, compute + * the averages of the Hessians of the function at the quadrature points, + * and return a vector with the correct type deduced by the interface + * Extractor you passed as the + * @p variable argument. + * + * Before you can call this method, you need to call the + * extract_local_dof_values() method at least once, passing the same + * @p global_vector_name string, and the same type for the variable @p dummy. + * + * If you have not previously called the extract_local_dof_values() method, + * this function will throw an exception. + * + * For this function to work properly, the underlying FEInterfaceValues + * object for which you called one of the reinit() functions, must have + * computed the information you are requesting. To do so, the + * update_hessians flag must be an element of the list of UpdateFlags that + * you passed to the constructor of this object. See "The interplay of + * UpdateFlags, Mapping, and FiniteElement" in the documentation of the + * FEValues class for more information. + */ + template + const std::vector< + typename FEInterfaceViews::View:: + template solution_hessian_type> & + get_averages_of_hessians(const std::string &global_vector_name, + const Extractor & variable, + const Number dummy = Number(0)); + + /** @} */ // CurrentInterfaceAverageEvaluation private: /** @@ -1016,8 +1246,7 @@ namespace MeshWorker const std::string &global_vector_name, Number dummy) const { - const unsigned int n_dofs = - get_current_fe_values().get_fe().n_dofs_per_cell(); + const unsigned int n_dofs = local_dof_indices.size(); const std::string dofs_name = get_unique_dofs_name(global_vector_name, n_dofs, dummy); @@ -1330,6 +1559,297 @@ namespace MeshWorker return ret; } + + + template + template + const std::vector:: + template solution_value_type> & + ScratchData::get_jumps_in_values( + const std::string &global_vector_name, + const Extractor & variable, + const Number dummy) + { + const std::vector &independent_local_dofs = + get_local_dof_values(global_vector_name, dummy); + + const FEInterfaceValues &feiv = + get_current_interface_fe_values(); + + AssertDimension(independent_local_dofs.size(), + feiv.n_current_interface_dofs()); + + const unsigned int n_q_points = feiv.n_quadrature_points; + + const std::string name = get_unique_name( + global_vector_name, variable, "_jump_values_q", n_q_points, dummy); + + // Now build the return type + using RetType = + std::vector:: + template solution_value_type>; + + RetType &ret = + internal_data_storage.template get_or_add_object_with_name( + name, n_q_points); + + AssertDimension(ret.size(), n_q_points); + + feiv[variable].get_jump_in_function_values_from_local_dof_values( + independent_local_dofs, ret); + return ret; + } + + + + template + template + const std::vector:: + template solution_gradient_type> & + ScratchData::get_jumps_in_gradients( + const std::string &global_vector_name, + const Extractor & variable, + const Number dummy) + { + const std::vector &independent_local_dofs = + get_local_dof_values(global_vector_name, dummy); + + const FEInterfaceValues &feiv = + get_current_interface_fe_values(); + + AssertDimension(independent_local_dofs.size(), + feiv.n_current_interface_dofs()); + + const unsigned int n_q_points = feiv.n_quadrature_points; + + const std::string name = get_unique_name( + global_vector_name, variable, "_jump_gradients_q", n_q_points, dummy); + + // Now build the return type + using RetType = + std::vector:: + template solution_gradient_type>; + + RetType &ret = + internal_data_storage.template get_or_add_object_with_name( + name, n_q_points); + + AssertDimension(ret.size(), n_q_points); + + feiv[variable].get_jump_in_function_gradients_from_local_dof_values( + independent_local_dofs, ret); + return ret; + } + + + + template + template + const std::vector:: + template solution_hessian_type> & + ScratchData::get_jumps_in_hessians( + const std::string &global_vector_name, + const Extractor & variable, + const Number dummy) + { + const std::vector &independent_local_dofs = + get_local_dof_values(global_vector_name, dummy); + + const FEInterfaceValues &feiv = + get_current_interface_fe_values(); + + AssertDimension(independent_local_dofs.size(), + feiv.n_current_interface_dofs()); + + const unsigned int n_q_points = feiv.n_quadrature_points; + + const std::string name = get_unique_name( + global_vector_name, variable, "_jump_hessians_q", n_q_points, dummy); + + // Now build the return type + using RetType = + std::vector:: + template solution_hessian_type>; + + RetType &ret = + internal_data_storage.template get_or_add_object_with_name( + name, n_q_points); + + AssertDimension(ret.size(), n_q_points); + + feiv[variable].get_jump_in_function_hessians_from_local_dof_values( + independent_local_dofs, ret); + return ret; + } + + + + template + template + const std::vector:: + template solution_third_derivative_type> & + ScratchData::get_jumps_in_third_derivatives( + const std::string &global_vector_name, + const Extractor & variable, + const Number dummy) + { + const std::vector &independent_local_dofs = + get_local_dof_values(global_vector_name, dummy); + + const FEInterfaceValues &feiv = + get_current_interface_fe_values(); + + AssertDimension(independent_local_dofs.size(), + feiv.n_current_interface_dofs()); + + const unsigned int n_q_points = feiv.n_quadrature_points; + + const std::string name = get_unique_name(global_vector_name, + variable, + "_jump_third_derivatives_q", + n_q_points, + dummy); + + // Now build the return type + using RetType = + std::vector:: + template solution_third_derivative_type>; + + RetType &ret = + internal_data_storage.template get_or_add_object_with_name( + name, n_q_points); + + AssertDimension(ret.size(), n_q_points); + + feiv[variable].get_jump_in_function_third_derivatives_from_local_dof_values( + independent_local_dofs, ret); + return ret; + } + + + + template + template + const std::vector:: + template solution_value_type> & + ScratchData::get_averages_of_values( + const std::string &global_vector_name, + const Extractor & variable, + const Number dummy) + { + const std::vector &independent_local_dofs = + get_local_dof_values(global_vector_name, dummy); + + const FEInterfaceValues &feiv = + get_current_interface_fe_values(); + + AssertDimension(independent_local_dofs.size(), + feiv.n_current_interface_dofs()); + + const unsigned int n_q_points = feiv.n_quadrature_points; + + const std::string name = get_unique_name( + global_vector_name, variable, "_average_values_q", n_q_points, dummy); + + // Now build the return type + using RetType = + std::vector:: + template solution_value_type>; + + RetType &ret = + internal_data_storage.template get_or_add_object_with_name( + name, n_q_points); + + AssertDimension(ret.size(), n_q_points); + + feiv[variable].get_average_of_function_values_from_local_dof_values( + independent_local_dofs, ret); + return ret; + } + + + + template + template + const std::vector:: + template solution_gradient_type> & + ScratchData::get_averages_of_gradients( + const std::string &global_vector_name, + const Extractor & variable, + const Number dummy) + { + const std::vector &independent_local_dofs = + get_local_dof_values(global_vector_name, dummy); + + const FEInterfaceValues &feiv = + get_current_interface_fe_values(); + + AssertDimension(independent_local_dofs.size(), + feiv.n_current_interface_dofs()); + + const unsigned int n_q_points = feiv.n_quadrature_points; + + const std::string name = get_unique_name( + global_vector_name, variable, "_average_gradients_q", n_q_points, dummy); + + // Now build the return type + using RetType = + std::vector:: + template solution_gradient_type>; + + RetType &ret = + internal_data_storage.template get_or_add_object_with_name( + name, n_q_points); + + AssertDimension(ret.size(), n_q_points); + + feiv[variable].get_average_of_function_gradients_from_local_dof_values( + independent_local_dofs, ret); + return ret; + } + + + + template + template + const std::vector:: + template solution_hessian_type> & + ScratchData::get_averages_of_hessians( + const std::string &global_vector_name, + const Extractor & variable, + const Number dummy) + { + const std::vector &independent_local_dofs = + get_local_dof_values(global_vector_name, dummy); + + const FEInterfaceValues &feiv = + get_current_interface_fe_values(); + + AssertDimension(independent_local_dofs.size(), + feiv.n_current_interface_dofs()); + + const unsigned int n_q_points = feiv.n_quadrature_points; + + const std::string name = get_unique_name( + global_vector_name, variable, "_average_hessians_q", n_q_points, dummy); + + // Now build the return type + using RetType = + std::vector:: + template solution_hessian_type>; + + RetType &ret = + internal_data_storage.template get_or_add_object_with_name( + name, n_q_points); + + AssertDimension(ret.size(), n_q_points); + + feiv[variable].get_average_of_function_hessians_from_local_dof_values( + independent_local_dofs, ret); + return ret; + } + + #endif } // namespace MeshWorker diff --git a/source/meshworker/scratch_data.cc b/source/meshworker/scratch_data.cc index 5198ebbd63..b5cf1e7ed2 100644 --- a/source/meshworker/scratch_data.cc +++ b/source/meshworker/scratch_data.cc @@ -291,6 +291,18 @@ namespace MeshWorker + template + const FEInterfaceValues & + ScratchData::get_current_interface_fe_values() const + { + Assert(interface_fe_values != nullptr, + ExcMessage("You have to initialize the cache using one of the " + "reinit functions first!")); + return *interface_fe_values; + } + + + template const FEValuesBase & ScratchData::get_current_neighbor_fe_values() const -- 2.39.5