From: Luca Heltai Date: Tue, 9 Apr 2019 21:04:36 +0000 (+0200) Subject: ScratchData X-Git-Tag: v9.1.0-rc1~210^2~10 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=5c7445b1ca3c72c36553ab859af1e0be2876de5d;p=dealii.git ScratchData --- diff --git a/include/deal.II/meshworker/scratch_data.h b/include/deal.II/meshworker/scratch_data.h new file mode 100644 index 0000000000..c249357c08 --- /dev/null +++ b/include/deal.II/meshworker/scratch_data.h @@ -0,0 +1,1264 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2019 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. +// +// --------------------------------------------------------------------- + +#ifndef dealii_meshworker_scratch_data_h +#define dealii_meshworker_scratch_data_h + +#include + +#include + +#include + +#include +#include + +#include + +#include +#include +#include + +DEAL_II_NAMESPACE_OPEN + +namespace MeshWorker +{ + namespace internal + { + /** + * Alias used to deduce the OutputType when asking for values or gradients + * of shape functions multiplied with the given @p NumberType. + */ + template + using OutputType = typename FEValuesViews::View:: + template OutputType; + } // namespace internal + + /** + * A helper class to simplify the parallel assembly of linear and non-linear + * problems, and the evaluation of finite element fields. + * + * This class is a drop in ScratchData to use with the WorkStream::run() + * function and with the MeshWorker::mesh_loop function(). + * + * The ScratchData class has three main goals: + * - to create FEValues, FEFaceValues, and FESubfaceValues for the current + * cell and for its neighbor cell *on demand* (only if they are + * necessary for the algorithm provided by the user), and to provide a + * uniform interface to access the FEValues objects when assembling cell, + * face, or subface contributions + * - to store arbitrary data types (or references to arbitrary data types), + * that can be retrieved by name (for example, when the assembly needs to + * make reference to other objects such as previous time steps' solution + * vectors, previous nonlinear iteration vectors, geometry descriptions, + * etc.) in an object of type + * - to provide a reasonable interface for those use cases where the user may + * need temporary vectors of data at quadrature points, allowing the + * construction of these temporaries *on demand*, and easy access to + * values, gradients, etc., of already computed solution vectors. + * + * The methods in @ref CurrentCellMethods initialize on demand internal FEValues, + * FEFaceValues, and FESubfaceValues objects on the current cell, allowing the + * use of this class as a single substitute for three different objects used + * to integrate and query finite element values on cells, faces, and subfaces. + * + * Similarly, the methods in @ref NeighborCellMethods initialize on demand + * (different) internal FEValues, FEFaceValues, and FESubfaceValues, on + * neighbor cells, allowing the use of this class aslo as a single substitute + * for the additional three objects you would typically need to integrate on + * the neighbor cell, and on its faces and subfaces (for example, in + * discontinuous Galerkin methods). + * + * If you need to retrieve values or gradients of finite element solution + * vectors, on the cell, face, or subface that has just beeing initialized + * with one of the functions in @ref CurrentCellMethods, you can use the + * methods in @ref CurrentCellEvaluation. + * + * An example usage for this class is given by the following snippet of code: + * + * @code + * ScratchData scratch(fe, + * cell_quadrature, + * update_values | update_gradients); + * + * FEValuesExtractors::Vector velocity(0); + * FEValuesExtractors::Scalar pressure(dim); + * + * ... + * + * for(const auto &cell : dof_handler.active_cell_iterators()) + * { + * const auto &fe_values = scratch.reinit(cell); + * scratch.extract_local_dof_values("current_solution", + * current_solution); + * + * scratch.extract_local_dof_values("previous_solution", + * previous_solution); + * + * const auto &local_solution_values = + * scratch.get_local_dof_values("current_solution", + * current_solution); + * + * const auto &local_previous_solution_values = + * scratch.get_local_dof_values("previous_solution", + * previous_solution); + * + * const auto &previous_solution_values_velocity = + * scratch.get_values("previous_solution", velocity); + * const auto &previous_solution_values_pressure = + * scratch.get_values("previous_solution", pressure); + * + * const auto &solution_values_velocity = + * scratch.get_values("solution", velocity); + * const auto &solution_values_pressure = + * scratch.get_values("solution", pressure); + * + * const auto &solution_symmetric_gradients_velocity = + * scratch.get_symmetric_gradients("solution", velocity); + * + * // Do something with the above. + * } + * @endcode + * + * The order in which you call functions of this class matters: if you call + * the ScratchData::reinit() function that takes an active cell iterator, + * then subsequent calls to methods that internally need an FEValuesBase + * object will use the internal FEValues object initialized with the given + * cell to perform their calculations. On the other hand, if you have called + * the ScratchData::reinit() method that also takes a face index, all + * subsequent calls to methods that need an FEValuesBase object, will use an + * internally stored FEFaceValues object, initialized with the cell and face + * index passed to the ScratchData::reinit() function. The same applies for + * the ScratchData::reinit() method that takes three arguments: the cell, the + * face index, and the subface index. + * + * The user code should be structured without interleaving work on cells and + * work on faces. + * + * Consider, for example, the following snippet of code: + * + * @code + * ScratchData scratch(...); + * FEValuesExtractors::Scalar temperature(0); + * + * for(const auto &cell : dof_handler.active_cell_iterators()) + * { + * const auto &fe_values = scratch.reinit(cell); + * const auto &local_dof_values = + * scratch.extract_local_dof_values("solution", solution); + * + * // This will contain all values of the temperature on the cell + * // quadrature points + * const auto &solution_values_cell = + * scratch.get_values("solution", temperature); + * + * // Do something with values on the cell + * ... + * + * // Now start working on the faces + * for(unsigned int f=0; f::faces_per_cell; ++f) + * { + * const auto &fe_face_values = scratch.reinit(cell, f); + * + * // Now we'll have access to the values of the temperature on the faces + * const auto &solution_values_face = + * scratch.get_values("solution", temperature); + * + * // Notice how the function call is the same, but the result depends + * // on what was the last reinit() function that was called. In this + * // case, we called reinit(cell, f), triggering internal work on + * // faces. + * } + * + * // At this point, we would like to go back and work on cells, + * // for example querying the values of the gradients: + * const auto &solution_gradients = + * scratch.get_gradients("solution", temperature); + * + * // This assertion would be triggered in debug mode + * AssertDimension(solution_gradients.size(), quadrature_cell.size()); + * + * // However, with the above call the content of solution_gradients is the + * // gradient of the temperature on the quadrature points of the last face + * // visited in the previous loop. + * // If you really want to have the values of the gradients on the cell, + * // at this point your only option is to call + * scratch.reinit(cell); + * + * // (again) before querying for the gradients: + * const auto &solution_gradients_cell = + * scratch.get_gradients("solution", temperature); + * + * // The call to reinit(cell), however, is an expensive one. You + * // should make sure you only call it once per cell by grouping together + * // all queries that concern cells in the same place (right after you + * // call the reinit(cell) method). + * // A similar argument holds for all calls on each face and subface. + * } + * @endcode + * + * When using this class, please cite + * @code{.bib} + * @article{SartoriGiulianiBardelloni-2018-a, + * Author = {Sartori, Alberto and Giuliani, Nicola and + * Bardelloni, Mauro and Heltai, Luca}, + * Journal = {SoftwareX}, + * Pages = {318--327}, + * Title = {{deal2lkit: A toolkit library for high performance + * programming in deal.II}}, + * Doi = {10.1016/j.softx.2018.09.004}, + * Volume = {7}, + * Year = {2018}} + * @endcode + * + * @author Luca Heltai, 2019. + */ + template + class ScratchData + { + public: + /** + * Create an empty ScratchData object. A SmartPointer pointing to + * @p mapping and @p fe is stored internally. Make sure they live longer + * than this class instance. + * + * The constructor does not initialize any of the internal FEValues objects. + * These are initialized the first time one of the reinit() functions is + * called, using the arguments passed here. + * + * @param mapping The mapping to use in the internal FEValues objects + * @param fe The finite element + * @param quadrature The cell quadrature + * @param update_flags UpdateFlags for the current cell FEValues and + * neighbor cell FEValues + * @param face_quadrature Face quadrature, used for FEFaceFvalues and + * FESubfaceValues for both the current cell and the neighbor cell + * @param face_update_flags UpdateFlags used for FEFaceFvalues and + * FESubfaceValues for both the current cell and the neighbor cell + */ + ScratchData( + const Mapping & mapping, + const FiniteElement &fe, + const Quadrature & quadrature, + const UpdateFlags & update_flags, + const Quadrature &face_quadrature = Quadrature(), + const UpdateFlags & face_update_flags = update_default); + + /** + * Similar to the other constructor, but this one allows to specify + * different flags for neighbor cells and faces. + * + * @param mapping The mapping to use in the internal FEValues objects + * @param fe The finite element + * @param quadrature The cell quadrature + * @param update_flags UpdateFlags for the current cell FEValues + * @param neighbor_update_flags UpdateFlags for the neighbor cell FEValues + * @param face_quadrature Face quadrature, used for FEFaceFvalues and + * FESubfaceValues for both the current cell and the neighbor cell + * @param face_update_flags UpdateFlags used for FEFaceFvalues and + * FESubfaceValues for the current cell + * @param neighbor_face_update_flags UpdateFlags used for FEFaceFvalues and + * FESubfaceValues for the neighbor cell + */ + ScratchData( + const Mapping & mapping, + const FiniteElement &fe, + const Quadrature & quadrature, + const UpdateFlags & update_flags, + const UpdateFlags & neighbor_update_flags, + const Quadrature &face_quadrature = Quadrature(), + const UpdateFlags & face_update_flags = update_default, + const UpdateFlags & neighbor_face_update_flags = update_default); + + /** + * Same as the other constructor, using the default MappingQ1. + * + * @param mapping The mapping to use in the internal FEValues objects + * @param fe The finite element + * @param quadrature The cell quadrature + * @param update_flags UpdateFlags for the current cell FEValues + * @param neighbor_update_flags UpdateFlags for the neighbor cell FEValues + * @param face_quadrature Face quadrature, used for FEFaceFvalues and + * FESubfaceValues for both the current cell and the neighbor cell + * @param face_update_flags UpdateFlags used for FEFaceFvalues and + * FESubfaceValues for the current cell + * @param neighbor_face_update_flags UpdateFlags used for FEFaceFvalues and + * FESubfaceValues for the neighbor cell + */ + ScratchData( + const FiniteElement &fe, + const Quadrature & quadrature, + const UpdateFlags & update_flags, + const Quadrature &face_quadrature = Quadrature(), + const UpdateFlags & face_update_flags = update_default); + + /** + * Same as the other constructor, using the default MappingQ1. + * + * @param mapping The mapping to use in the internal FEValues objects + * @param fe The finite element + * @param quadrature The cell quadrature + * @param update_flags UpdateFlags for the current cell FEValues + * @param neighbor_update_flags UpdateFlags for the neighbor cell FEValues + * @param face_quadrature Face quadrature, used for FEFaceFvalues and + * FESubfaceValues for both the current cell and the neighbor cell + * @param face_update_flags UpdateFlags used for FEFaceFvalues and + * FESubfaceValues for the current cell + * @param neighbor_face_update_flags UpdateFlags used for FEFaceFvalues and + * FESubfaceValues for the neighbor cell + */ + ScratchData( + const FiniteElement &fe, + const Quadrature & quadrature, + const UpdateFlags & update_flags, + const UpdateFlags & neighbor_update_flags, + const Quadrature &face_quadrature = Quadrature(), + const UpdateFlags & face_update_flags = update_default, + const UpdateFlags & neighbor_face_update_flags = update_default); + + /** + * Deep copy constructor. FEValues objects are not copied. + */ + ScratchData(const ScratchData &scratch); + + /** + * \defgroup CurrentCellMethos Methods to work on current cell + */ + /**@{*/ + + /** + * Initialize the internal FEValues with the given @p cell, and return + * a reference to it. + * + * After calling this function, get_current_fe_values() will return the + * same object of this method, as an FEValuesBase reference. + */ + const FEValues & + reinit( + const typename DoFHandler::active_cell_iterator &cell); + + /** + * Initialize the internal FEFaceValues to use the given @p face_no on the given + * @p cell, and return a reference to it. + * + * After calling this function, get_current_fe_values() will return the + * same object of this method, as an FEValuesBase reference. + */ + const FEFaceValues & + reinit(const typename DoFHandler::active_cell_iterator &cell, + const unsigned int face_no); + + /** + * Initialize the internal FESubfaceValues to use the given @p subface_no, + * on @p face_no, on the given @p cell, and return a reference to it. + * + * After calling this function, get_current_fe_values() will return the + * same object of this method, as an FEValuesBase reference. + * + * If @p subface_no is numbers::invalid_unsigned_int, the reinit() function + * that takes only the @p cell and the @p face_no is called. + */ + const FEFaceValuesBase & + reinit(const typename DoFHandler::active_cell_iterator &cell, + const unsigned int face_no, + const unsigned int subface_no); + + /** + * Get the currently initialized FEValues. + * + * This function will return the internal FEValues if the + * reinit(cell) function was called last. If the reinit(cell, face_no) + * function was called, then this function returns the internal + * FEFaceValues, and if the reinit(cell, face_no, subface_no) function was + * called (with a valid @p subface_no argument), it returns the internal + * FESubfaceValues object. + */ + const FEValuesBase & + get_current_fe_values() const; + + /** + * Return the quadrature points of the internal FEValues object. + */ + const std::vector> & + get_quadrature_points() const; + + /** + * Return the JxW values of the internal FEValues object. + */ + const std::vector & + get_JxW_values() const; + + /** + * Return the last computed normal vectors. + */ + const std::vector> & + get_normal_vectors() const; + + /** + * Return the local dof indices for the cell passed the last time the + * reinit() function was called. + */ + const std::vector & + get_local_dof_indices() const; + + /** @} */ // CurrentCellMethods + + /** + * @defgroup NeighborCellMethods Methods to work on neighbor cell + */ + /** @{ */ + + /** + * Initialize the internal neighbor FEValues to use the given @p cell, and + * return a reference to it. + * + * After calling this function, get_current_fe_values() will return the + * same object of this method, as an FEValuesBase reference. + */ + const FEValues & + reinit_neighbor( + const typename DoFHandler::active_cell_iterator &cell); + + /** + * Initialize the internal FEFaceValues to use the given @p face_no on the + * given @p cell, and return a reference to it. + * + * After calling this function, get_current_fe_values() will return the + * same object of this method, as an FEValuesBase reference. + */ + const FEFaceValues & + reinit_neighbor( + const typename DoFHandler::active_cell_iterator &cell, + const unsigned int face_no); + + /** + * Initialize the internal FESubfaceValues to use the given @p subface_no, + * on @p face_no, on the given @p cell, and return a reference to it. + * + * After calling this function, get_current_fe_values() will return the + * same object of this method, as an FEValuesBase reference. + * + * If @p subface_no is numbers::invalid_unsigned_int, the reinit() function + * that takes only the @p cell and the @p face_no is called. + */ + const FEFaceValuesBase & + reinit_neighbor( + const typename DoFHandler::active_cell_iterator &cell, + const unsigned int face_no, + const unsigned int subface_no); + + /** + * Get the currently initialized neighbor FEValues. + * + * This function will return the neighbor FEValues if the + * reinit_neighbor(cell) function was called last. If the + * reinit_neighbor(cell, face_no) function was called, then this function + * returns the internal neighbor FEFaceValues, and if the + * reinit_neighbor(cell, face_no, subface_no) function was + * called (with a valid @p subface_no argument), it returns the internal neighbor + * FESubfaceValues object. + */ + const FEValuesBase & + get_current_neighbor_fe_values() const; + + /** + * Return the JxW values of the neighbor FEValues object. + */ + const std::vector & + get_neighbor_JxW_values() const; + + /** + * Return the last computed normal vectors on the neighbor. + */ + const std::vector> & + get_neighbor_normal_vectors(); + + /** + * Return the local dof indices of the neighbor passed the last time the + * reinit_neighbor() function was called. + */ + const std::vector & + get_neighbor_dof_indices() const; + + /** @} */ // NeighborCellMethods + + /** + * Return a GeneralDataStorage object that can be used to store store any + * amount of data, of any type, which is then made accessible by an + * identifier string. + */ + GeneralDataStorage & + get_general_data_storage(); + + /** @defgroup CurrentCellEvaluation Evaluation of finite element fields + * and their derivatives on the current cell + * @{ + */ + + /** + * Extract the local dof values associated with the internally initialized + * cell. + * + * Before you call this function, you have to make sure you have previously + * called one of the reinit() functions. + * + * At every call of this function, a new vector of dof values is generated + * and stored internally, unless a previous vector with the same name is + * found. If this is the case, the content of that vector is overwritten. + * + * If you give a unique @p global_vector_name, then for each cell you are + * guaranteed to get a unique vector of independent dofs of the same type + * as the dummy variable. If you use an automatic differentiation number + * type (like Sacado::Fad::DFad, + * Sacado::Fad::DFad>, etc.) this method will + * also initialize the independent variables internally, allowing you + * to perform automatic differentiation. + * + * You can access the extracted local dof values by calling the method + * get_local_dof_values() with the same @p global_vector_name argument + * you passed here. + * + * Notice that using this initialization strategy renders the use of this + * ScratchData object incompatible with the AD helper classes (since they + * do their own data management). In particular, it is necessary for the + * user to manage all of the AD data themselves (both before and after this + * call). + */ + template + void + extract_local_dof_values(const std::string &global_vector_name, + const VectorType & input_vector, + const Number dummy = Number(0)); + + /** + * After calling extract_local_dof_values(), you can retrieve the stored + * information through this method. + * + * Both the argument @p global_vector_name and the type of the @p dummy + * variable should match the ones you passed to the + * extract_local_dof_values() function. + */ + template + const std::vector & + get_local_dof_values(const std::string &global_vector_name, + Number dummy = Number(0)) const; + + /** + * For the solution vector identified by @p global_vector_name, compute + * the values 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. + * + * 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 FEValues, + * FEFaceValues, or FESubfaceValues 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 internal::OutputType:: + value_type> & + get_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 gradients 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. + * + * 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 FEValues, + * FEFaceValues, or FESubfaceValues 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 internal::OutputType:: + gradient_type> & + get_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 symmetric gradients 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. + * + * 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 FEValues, + * FEFaceValues, or FESubfaceValues 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 internal::OutputType:: + symmetric_gradient_type> & + get_symmetric_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 divergences 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. + * + * 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 FEValues, + * FEFaceValues, or FESubfaceValues 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 internal::OutputType:: + divergence_type> & + get_divergences(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 curls 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. + * + * 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 FEValues, + * FEFaceValues, or FESubfaceValues 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::curl_type> + & + get_curls(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 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. + * + * 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 FEValues, + * FEFaceValues, or FESubfaceValues 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 internal::OutputType:: + hessian_type> & + get_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 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. + * + * 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 FEValues, + * FEFaceValues, or FESubfaceValues 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 internal::OutputType:: + third_derivative_type> & + get_third_derivatives(const std::string &global_vector_name, + const Extractor & variable, + const Number dummy = Number(0)); + + /** @} */ // CurrentCellEvaluation + + /** + * Return a reference to the used mapping. + */ + const Mapping & + get_mapping() const; + + private: + /** + * Construct a unique name to store vectors of values, gradients, + * divergences, etc., in the internal GeneralDataStorage object. + */ + template + std::string + get_unique_name(const std::string &global_vector_name, + const Extractor & variable, + const std::string &object_type, + const unsigned int size, + const Number & exemplar_number) const; + + /** + * Construct a unique name to store local dof values. + */ + template + std::string + get_unique_dofs_name(const std::string &global_vector_name, + const unsigned int size, + const Number & exemplar_number) const; + + /** + * The mapping used by the internal FEValues. Make sure it lives + * longer than this class. + */ + SmartPointer> mapping; + + /** + * The finite element used by the internal FEValues. Make sure it lives + * longer than this class. + */ + SmartPointer> fe; + + /** + * Quadrature formula used to integrate on the current cell, and on its + * neighbor. + */ + Quadrature cell_quadrature; + + /** + * Quadrature formula used to integrate on faces, subfaces, and neighbor + * faces and subfaces. + */ + Quadrature face_quadrature; + + /** + * UpdateFlags to use when initializing the cell FEValues object. + */ + UpdateFlags cell_update_flags; + + /** + * UpdateFlags to use when initializing the neighbor cell FEValues objects. + */ + UpdateFlags neighbor_cell_update_flags; + + /** + * UpdateFlags to use when initializing FEFaceValues and FESubfaceValues + * objects. + */ + UpdateFlags face_update_flags; + + /** + * UpdateFlags to use when initializing neighbor FEFaceValues and + * FESubfaceValues objects. + */ + UpdateFlags neighbor_face_update_flags; + + /** + * Finite element values on the current cell. + */ + std::unique_ptr> fe_values; + + /** + * Finite element values on the current face. + */ + std::unique_ptr> fe_face_values; + + /** + * Finite element values on the current subface. + */ + std::unique_ptr> fe_subface_values; + + /** + * Finite element values on the neighbor cell. + */ + std::unique_ptr> neighbor_fe_values; + + /** + * Finite element values on the neighbor face. + */ + std::unique_ptr> neighbor_fe_face_values; + + /** + * Finite element values on the neighbor subface. + */ + std::unique_ptr> neighbor_fe_subface_values; + + /** + * Dof indices on the current cell. + */ + std::vector local_dof_indices; + + /** + * Dof indices on the neighbor cell. + */ + std::vector neighbor_dof_indices; + + /** + * User data storage. + */ + GeneralDataStorage user_data_storage; + + /** + * Internal data storage. + */ + GeneralDataStorage internal_data_storage; + + /** + * A pointer to the last used FEValues/FEFaceValues, or FESubfaceValues + * object on this cell. + */ + SmartPointer> current_fe_values; + + /** + * A pointer to the last used FEValues/FEFaceValues, or FESubfaceValues + * object on the neighbor cell. + */ + SmartPointer> current_neighbor_fe_values; + }; + +#ifndef DOXYGEN + template + template + std::string + ScratchData::get_unique_name( + const std::string &global_vector_name, + const Extractor & variable, + const std::string &object_type, + const unsigned int size, + const Number & exemplar_number) const + { + return global_vector_name + "_" + variable.get_name() + "_" + object_type + + "_" + Utilities::int_to_string(size) + "_" + + Utilities::type_to_string(exemplar_number); + } + + + + template + template + std::string + ScratchData::get_unique_dofs_name( + const std::string &global_vector_name, + const unsigned int size, + const Number & exemplar_number) const + { + return global_vector_name + "_independent_local_dofs_" + + Utilities::int_to_string(size) + "_" + + Utilities::type_to_string(exemplar_number); + } + + + + template + template + void + ScratchData::extract_local_dof_values( + const std::string &global_vector_name, + const VectorType & input_vector, + const Number dummy) + { + const unsigned int n_dofs = get_current_fe_values().get_fe().dofs_per_cell; + + const std::string name = + get_unique_dofs_name(global_vector_name, n_dofs, dummy); + + auto &independent_local_dofs = + internal_data_storage + .template get_or_add_object_with_name>(name, + n_dofs); + + AssertDimension(independent_local_dofs.size(), n_dofs); + + if (Differentiation::AD::is_tapeless_ad_number::value == true) + for (unsigned int i = 0; i < n_dofs; ++i) + Differentiation::AD::internal::Marking::independent_variable( + input_vector(local_dof_indices[i]), + i, + n_dofs, + independent_local_dofs[i]); + else + for (unsigned int i = 0; i < n_dofs; ++i) + independent_local_dofs[i] = input_vector(local_dof_indices[i]); + } + + + + template + template + const std::vector & + ScratchData::get_local_dof_values( + const std::string &global_vector_name, + Number dummy) const + { + const unsigned int n_dofs = get_current_fe_values().get_fe().dofs_per_cell; + + const std::string dofs_name = + get_unique_dofs_name(global_vector_name, n_dofs, dummy); + + Assert( + internal_data_storage.stores_object_with_name(dofs_name), + ExcMessage( + "You did not call yet extract_local_dof_values with the right types!")); + + return internal_data_storage + .template get_object_with_name>(dofs_name); + } + + + + template + template + const std::vector< + typename internal::OutputType::value_type> + & + ScratchData::get_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 FEValuesBase &fev = get_current_fe_values(); + + const unsigned int n_q_points = fev.n_quadrature_points; + + const std::string name = get_unique_name( + global_vector_name, variable, "_values_q", n_q_points, dummy); + + // Now build the return type + using RetType = + std::vector::value_type>; + + RetType &ret = + internal_data_storage.template get_or_add_object_with_name( + name, n_q_points); + + AssertDimension(ret.size(), n_q_points); + + fev[variable].get_function_values_from_local_dof_values( + independent_local_dofs, ret); + return ret; + } + + + + template + template + const std::vector< + typename internal::OutputType:: + gradient_type> & + ScratchData::get_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 FEValuesBase &fev = get_current_fe_values(); + + const unsigned int n_q_points = fev.n_quadrature_points; + + const std::string name = get_unique_name( + global_vector_name, variable, "_gradients_q", n_q_points, dummy); + + // Now build the return type + using RetType = std::vector< + typename internal::OutputType:: + gradient_type>; + + RetType &ret = + internal_data_storage.template get_or_add_object_with_name( + name, n_q_points); + + AssertDimension(ret.size(), n_q_points); + + fev[variable].get_function_gradients_from_local_dof_values( + independent_local_dofs, ret); + return ret; + } + + + + template + template + const std::vector< + typename internal::OutputType:: + hessian_type> & + ScratchData::get_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 FEValuesBase &fev = get_current_fe_values(); + + const unsigned int n_q_points = fev.n_quadrature_points; + + const std::string name = get_unique_name( + global_vector_name, variable, "_hessians_q", n_q_points, dummy); + + // Now build the return type + using RetType = + std::vector::hessian_type>; + + RetType &ret = + internal_data_storage.template get_or_add_object_with_name( + name, n_q_points); + + + AssertDimension(ret.size(), n_q_points); + + fev[variable].get_function_hessians_from_local_dof_values( + independent_local_dofs, ret); + return ret; + } + + + + template + template + const std::vector< + typename internal::OutputType:: + third_derivative_type> & + ScratchData::get_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 FEValuesBase &fev = get_current_fe_values(); + + const unsigned int n_q_points = fev.n_quadrature_points; + + const std::string name = get_unique_name( + global_vector_name, variable, "_third_derivatives_q", n_q_points, dummy); + + // Now build the return type + using RetType = std::vector< + typename internal::OutputType:: + 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); + + fev[variable].get_function_third_derivatives_from_local_dof_values( + independent_local_dofs, ret); + return ret; + } + + + + template + template + const std::vector< + typename internal::OutputType:: + symmetric_gradient_type> & + ScratchData::get_symmetric_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 FEValuesBase &fev = get_current_fe_values(); + + const unsigned int n_q_points = fev.n_quadrature_points; + + const std::string name = get_unique_name( + global_vector_name, variable, "_symmetric_gradient_q", n_q_points, dummy); + + + // Now build the return type + using RetType = std::vector< + typename internal::OutputType:: + symmetric_gradient_type>; + + RetType &ret = + internal_data_storage.template get_or_add_object_with_name( + name, n_q_points); + + + AssertDimension(ret.size(), n_q_points); + + fev[variable].get_function_symmetric_gradients_from_local_dof_values( + independent_local_dofs, ret); + return ret; + } + + + template + template + const std::vector< + typename internal::OutputType:: + divergence_type> & + ScratchData::get_divergences( + 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 FEValuesBase &fev = get_current_fe_values(); + + const unsigned int n_q_points = fev.n_quadrature_points; + + const std::string name = get_unique_name( + global_vector_name, variable, "_divergence_q", n_q_points, dummy); + + // Now build the return type + using RetType = std::vector< + typename internal::OutputType:: + divergence_type>; + + RetType &ret = + internal_data_storage.template get_or_add_object_with_name( + name, n_q_points); + + + AssertDimension(ret.size(), n_q_points); + + fev[variable].get_function_divergences_from_local_dof_values( + independent_local_dofs, ret); + return ret; + } + + + + template + template + const std::vector< + typename internal::OutputType::curl_type> + & + ScratchData::get_curls(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 FEValuesBase &fev = get_current_fe_values(); + + const unsigned int n_q_points = fev.n_quadrature_points; + + const std::string name = get_unique_name( + global_vector_name, variable, "_curl_q", n_q_points, dummy); + + // Now build the return type + using RetType = + std::vector::curl_type>; + + RetType &ret = + internal_data_storage.template get_or_add_object_with_name( + name, n_q_points); + + AssertDimension(ret.size(), n_q_points); + + fev[variable].get_function_curls_from_local_dof_values( + independent_local_dofs, ret); + return ret; + } + +#endif + +} // namespace MeshWorker + +DEAL_II_NAMESPACE_CLOSE + +#endif diff --git a/source/meshworker/CMakeLists.txt b/source/meshworker/CMakeLists.txt index 557829dbad..f939da4894 100644 --- a/source/meshworker/CMakeLists.txt +++ b/source/meshworker/CMakeLists.txt @@ -19,11 +19,13 @@ SET(_src mesh_worker.cc mesh_worker_info.cc mesh_worker_vector_selector.cc + scratch_data.cc ) SET(_inst mesh_worker_info.inst.in mesh_worker_vector_selector.inst.in + scratch_data.inst.in ) FILE(GLOB _header diff --git a/source/meshworker/scratch_data.cc b/source/meshworker/scratch_data.cc new file mode 100644 index 0000000000..8e6e8708c8 --- /dev/null +++ b/source/meshworker/scratch_data.cc @@ -0,0 +1,356 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2019 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. +// +// --------------------------------------------------------------------- + +#include + +DEAL_II_NAMESPACE_OPEN + +namespace MeshWorker +{ + template + ScratchData::ScratchData( + const Mapping & mapping, + const FiniteElement &fe, + const Quadrature & quadrature, + const UpdateFlags & update_flags, + const Quadrature & face_quadrature, + const UpdateFlags & face_update_flags) + : mapping(&mapping) + , fe(&fe) + , cell_quadrature(quadrature) + , face_quadrature(face_quadrature) + , cell_update_flags(update_flags) + , neighbor_cell_update_flags(update_flags) + , face_update_flags(face_update_flags) + , neighbor_face_update_flags(face_update_flags) + , local_dof_indices(fe.dofs_per_cell) + , neighbor_dof_indices(fe.dofs_per_cell) + {} + + + + template + ScratchData::ScratchData( + const Mapping & mapping, + const FiniteElement &fe, + const Quadrature & quadrature, + const UpdateFlags & update_flags, + const UpdateFlags & neighbor_update_flags, + const Quadrature & face_quadrature, + const UpdateFlags & face_update_flags, + const UpdateFlags & neighbor_face_update_flags) + : mapping(&mapping) + , fe(&fe) + , cell_quadrature(quadrature) + , face_quadrature(face_quadrature) + , cell_update_flags(update_flags) + , neighbor_cell_update_flags(neighbor_update_flags) + , face_update_flags(face_update_flags) + , neighbor_face_update_flags(neighbor_face_update_flags) + , local_dof_indices(fe.dofs_per_cell) + , neighbor_dof_indices(fe.dofs_per_cell) + {} + + + + template + ScratchData::ScratchData( + const FiniteElement &fe, + const Quadrature & quadrature, + const UpdateFlags & update_flags, + const Quadrature & face_quadrature, + const UpdateFlags & face_update_flags) + : ScratchData(StaticMappingQ1::mapping, + fe, + quadrature, + update_flags, + face_quadrature, + face_update_flags) + {} + + + + template + ScratchData::ScratchData( + const FiniteElement &fe, + const Quadrature & quadrature, + const UpdateFlags & update_flags, + const UpdateFlags & neighbor_update_flags, + const Quadrature & face_quadrature, + const UpdateFlags & face_update_flags, + const UpdateFlags & neighbor_face_update_flags) + : ScratchData(StaticMappingQ1::mapping, + fe, + quadrature, + update_flags, + neighbor_update_flags, + face_quadrature, + face_update_flags, + neighbor_face_update_flags) + {} + + + + template + ScratchData::ScratchData( + const ScratchData &scratch) + : mapping(scratch.mapping) + , fe(scratch.fe) + , cell_quadrature(scratch.cell_quadrature) + , face_quadrature(scratch.face_quadrature) + , cell_update_flags(scratch.cell_update_flags) + , neighbor_cell_update_flags(scratch.neighbor_cell_update_flags) + , face_update_flags(scratch.face_update_flags) + , neighbor_face_update_flags(scratch.neighbor_face_update_flags) + , local_dof_indices(scratch.local_dof_indices) + , neighbor_dof_indices(scratch.neighbor_dof_indices) + , user_data_storage(scratch.user_data_storage) + , internal_data_storage(scratch.internal_data_storage) + {} + + + + template + const FEValues & + ScratchData::reinit( + const typename DoFHandler::active_cell_iterator &cell) + { + if (!fe_values) + fe_values = std::make_unique>(*mapping, + *fe, + cell_quadrature, + cell_update_flags); + + fe_values->reinit(cell); + cell->get_dof_indices(local_dof_indices); + current_fe_values = fe_values.get(); + return *fe_values; + } + + + + template + const FEFaceValues & + ScratchData::reinit( + const typename DoFHandler::active_cell_iterator &cell, + const unsigned int face_no) + { + if (!fe_face_values) + fe_face_values = std::make_unique>( + *mapping, *fe, face_quadrature, face_update_flags); + + fe_face_values->reinit(cell, face_no); + cell->get_dof_indices(local_dof_indices); + current_fe_values = fe_face_values.get(); + return *fe_face_values; + } + + + + template + const FEFaceValuesBase & + ScratchData::reinit( + const typename DoFHandler::active_cell_iterator &cell, + const unsigned int face_no, + const unsigned int subface_no) + { + if (subface_no != numbers::invalid_unsigned_int) + { + if (!fe_subface_values) + fe_subface_values = std::make_unique>( + *mapping, *fe, face_quadrature, face_update_flags); + fe_subface_values->reinit(cell, face_no, subface_no); + cell->get_dof_indices(local_dof_indices); + + current_fe_values = fe_subface_values.get(); + return *fe_subface_values; + } + else + return reinit(cell, face_no); + } + + + + template + const FEValues & + ScratchData::reinit_neighbor( + const typename DoFHandler::active_cell_iterator &cell) + { + if (!neighbor_fe_values) + neighbor_fe_values = std::make_unique>( + *mapping, *fe, cell_quadrature, neighbor_cell_update_flags); + + neighbor_fe_values->reinit(cell); + cell->get_dof_indices(neighbor_dof_indices); + current_neighbor_fe_values = neighbor_fe_values.get(); + return *neighbor_fe_values; + } + + + + template + const FEFaceValues & + ScratchData::reinit_neighbor( + const typename DoFHandler::active_cell_iterator &cell, + const unsigned int face_no) + { + if (!neighbor_fe_face_values) + neighbor_fe_face_values = std::make_unique>( + *mapping, *fe, face_quadrature, neighbor_face_update_flags); + neighbor_fe_face_values->reinit(cell, face_no); + cell->get_dof_indices(neighbor_dof_indices); + current_neighbor_fe_values = neighbor_fe_face_values.get(); + return *neighbor_fe_face_values; + } + + + + template + const FEFaceValuesBase & + ScratchData::reinit_neighbor( + const typename DoFHandler::active_cell_iterator &cell, + const unsigned int face_no, + const unsigned int subface_no) + { + if (subface_no != numbers::invalid_unsigned_int) + { + if (!neighbor_fe_subface_values) + neighbor_fe_subface_values = + std::make_unique>( + *mapping, *fe, face_quadrature, neighbor_face_update_flags); + neighbor_fe_subface_values->reinit(cell, face_no, subface_no); + cell->get_dof_indices(neighbor_dof_indices); + current_neighbor_fe_values = neighbor_fe_subface_values.get(); + return *neighbor_fe_subface_values; + } + else + return reinit_neighbor(cell, face_no); + } + + + + template + const FEValuesBase & + ScratchData::get_current_fe_values() const + { + Assert(current_fe_values != nullptr, + ExcMessage("You have to initialize the cache using one of the " + "reinit functions first!")); + return *current_fe_values; + } + + + + template + const FEValuesBase & + ScratchData::get_current_neighbor_fe_values() const + { + Assert(current_neighbor_fe_values != nullptr, + ExcMessage("You have to initialize the cache using one of the " + "reinit functions first!")); + return *current_neighbor_fe_values; + } + + + + template + const std::vector> & + ScratchData::get_quadrature_points() const + { + return get_current_fe_values().get_quadrature_points(); + } + + + + template + const std::vector & + ScratchData::get_JxW_values() const + { + return get_current_fe_values().get_JxW_values(); + } + + + + template + const std::vector & + ScratchData::get_neighbor_JxW_values() const + { + return get_current_neighbor_fe_values().get_JxW_values(); + } + + + + template + const std::vector> & + ScratchData::get_normal_vectors() const + { + return get_current_fe_values().get_normal_vectors(); + } + + + + template + const std::vector> & + ScratchData::get_neighbor_normal_vectors() + { + return get_current_neighbor_fe_values().get_normal_vectors(); + } + + + + template + const std::vector & + ScratchData::get_local_dof_indices() const + { + return local_dof_indices; + } + + + + template + const std::vector & + ScratchData::get_neighbor_dof_indices() const + { + return neighbor_dof_indices; + } + + + + template + GeneralDataStorage & + ScratchData::get_general_data_storage() + { + return user_data_storage; + } + + + + template + const Mapping & + ScratchData::get_mapping() const + { + return *mapping; + } + +} // namespace MeshWorker +DEAL_II_NAMESPACE_CLOSE + +// Explicit instantiations +DEAL_II_NAMESPACE_OPEN +namespace MeshWorker +{ +#include "scratch_data.inst" +} +DEAL_II_NAMESPACE_CLOSE \ No newline at end of file diff --git a/source/meshworker/scratch_data.inst.in b/source/meshworker/scratch_data.inst.in new file mode 100644 index 0000000000..96eb79e413 --- /dev/null +++ b/source/meshworker/scratch_data.inst.in @@ -0,0 +1,23 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2019 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. +// +// --------------------------------------------------------------------- + + +for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS) + { +#if deal_II_dimension <= deal_II_space_dimension + template class ScratchData; + +#endif + }