From f0eb0721a132db569935a50bffc06d0947e6b875 Mon Sep 17 00:00:00 2001 From: Denis Davydov Date: Fri, 22 Apr 2016 17:45:05 +0200 Subject: [PATCH] add classes to manage quadrature point data, including transfering between cells --- doc/news/changes.h | 10 + include/deal.II/base/quadrature_point_data.h | 867 ++++++++++++++++++ tests/base/cell_data_storage_01.cc | 171 ++++ tests/base/cell_data_storage_01.output | 2 + tests/base/cell_data_storage_02.cc | 152 +++ tests/base/cell_data_storage_02.output | 2 + tests/base/quadrature_point_data.cc | 219 +++++ ...e_point_data.with_p4est=on.mpirun=3.output | 2 + tests/base/quadrature_point_data_02.cc | 233 +++++ ...oint_data_02.with_p4est=on.mpirun=3.output | 2 + tests/base/quadrature_point_data_03.cc | 238 +++++ ...oint_data_03.with_p4est=on.mpirun=3.output | 2 + 12 files changed, 1900 insertions(+) create mode 100644 include/deal.II/base/quadrature_point_data.h create mode 100644 tests/base/cell_data_storage_01.cc create mode 100644 tests/base/cell_data_storage_01.output create mode 100644 tests/base/cell_data_storage_02.cc create mode 100644 tests/base/cell_data_storage_02.output create mode 100644 tests/base/quadrature_point_data.cc create mode 100644 tests/base/quadrature_point_data.with_p4est=on.mpirun=3.output create mode 100644 tests/base/quadrature_point_data_02.cc create mode 100644 tests/base/quadrature_point_data_02.with_p4est=on.mpirun=3.output create mode 100644 tests/base/quadrature_point_data_03.cc create mode 100644 tests/base/quadrature_point_data_03.with_p4est=on.mpirun=3.output diff --git a/doc/news/changes.h b/doc/news/changes.h index 75d4d611d6..adbfeb3c43 100644 --- a/doc/news/changes.h +++ b/doc/news/changes.h @@ -120,6 +120,16 @@ inconvenience this causes.

General

    +
  1. New: Add a collection of classes to manage user's quadrature point data: + CellDataStorage, TransferableQuadraturePointData and + parallel::distributed::ContinuousQuadratureDataTransfer. + The implementation of CellDataStorage is flexible to support different types of + data object at different cells. parallel::distributed::ContinuousQuadratureDataTransfer + provides a convenient interface to transfer quadrature point data between cells + of parallel::distributed::Triangulation. +
    + (Denis Davydov, Jean-Paul Pelteret, 2016/04/30) +
  2. New: Added an interface to the GNU Scientific Library. Also introduce a cubic spline interpolation function Functions::CSpline. diff --git a/include/deal.II/base/quadrature_point_data.h b/include/deal.II/base/quadrature_point_data.h new file mode 100644 index 0000000000..251dde8e30 --- /dev/null +++ b/include/deal.II/base/quadrature_point_data.h @@ -0,0 +1,867 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2016 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + +#ifndef dealii__quadrature_point_data_h +#define dealii__quadrature_point_data_h + +#ifdef DEAL_II_WITH_CXX11 + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include + +DEAL_II_NAMESPACE_OPEN + +/*!@addtogroup Quadrature */ +/*@{*/ + +/** + * A class for storing at each cell represented by iterators of type @p CellIteratorType + * a vector of data @p DataType . + * The underlying structure and the initialize() method of this class are designed + * in such a way that one could use different child classes derived from the + * base DataType class to store data on a given cell. This implies the usage of pointers, + * in our case -- std::shared_ptr(). + * + * @note The data type stored on each cell can be different. + * However, within the cell this class stores a vector of objects of a single data type. + * For this reason, this class may not be sufficiently flexible when, for example, + * adopting a level-set approach to describe material behavior. + * + * @note This class is only available if deal.II is configured with C++11 support. + * + * @author Denis Davydov, Jean-Paul Pelteret, 2016 + */ +template +class CellDataStorage +{ +public: + /** + * A typedef for a vector of shared pointers with arbitrary data type @p T . + */ + template using CellVector = std::vector >; + + /** + * Default constructor. + */ + CellDataStorage(); + + /** + * Default destructor. + */ + ~CellDataStorage(); + + /** + * Initialize data on the @p cell to store @p number_of_data_points_per_cell of objects of type @p T . + * The typename @p T is possibly another class which is derived from the + * base @p DataType class. In order to initialize the vector of objects + * we have to assume that the class @p T has a default constructor. + * This function has to be called on every cell where data is to be stored. + * + * After the data is initialized, it can be modified using the get_data() + * function. + * + * @note Subsequent calls of this function with the same @p cell will not + * alter the objects associated with it. In order to remove the stored data, + * use the erase() function. + * + * @note It is possible to use different types @p T for different cells which + * may reflect, for example, different constitutive models of continuum mechanics + * in different parts of the domain. + * + * @pre The type @p T needs to either equal @p DataType, or be a class derived + * from @p DataType. + */ + template + void initialize(const CellIteratorType &cell, + const unsigned int number_of_data_points_per_cell); + + /** + * Same as above but for a range of iterators starting at @p cell_start + * until, but not including, @p cell_end for all locally owned cells, i.e. + * for which `cell->is_locally_owned()==true` . + */ + template + void initialize(const CellIteratorType &cell_start, + const CellIteratorType &cell_end, + const unsigned int number_of_data_points_per_cell); + + /** + * Removes data stored at the @p cell. Returns true if the data was removed. + * If no data is attached to the @p cell, this function will not do anything + * and returns false. + * + * @note This function will also check that there are no + * outstanding references to the data stored on this cell. That is to say, + * that the only references to the stored data are that made by this class. + */ + bool erase(const CellIteratorType &cell); + + /** + * Clear all the data stored in this object. + */ + void clear(); + + /** + * Get a vector of the data located at @p cell . + * A possible additional typename @p T is the class to which the base class + * DataType could be cast. Since @p DataType is stored as shared pointers, + * there is minimal overhead in returning a vector by value instead of by reference. + * This allows flexibility if class @p T is not the same as @p DataType on a + * cell-by-cell basis. + * + * @pre The type @p T needs to match the class provided to initialize() . + */ + template + std::vector > get_data(const CellIteratorType &cell); + + /** + * Get a vector of constant pointers to data located at @p cell . + * A possible additional typename @p T is the class to which the base class + * DataType could be cast. Since @p DataType is stored as shared pointers, + * there is minimal overhead in returning a vector by value instead of by reference. + * This allows flexibility if class @p T is not the same as @p DataType on a + * cell-by-cell basis. + * + * @pre The type @p T needs to match the class provided to initialize() . + */ + template + std::vector > get_data(const CellIteratorType &cell) const; + +private: + /** + * A map to store a vector of data on a cell. + */ + std::map > map; + +}; + + +/** + * An abstract class which specifies requirements for data on + * a single quadrature point to be transferable during refinement or repartitioning. + * + * This class provides a framework by which derived classes representing data at + * quadrature points can declare how many scalar values they store, and then + * implement functions that pack and unpack these scalars into arrays. + * These arrays are used to transfer data from quadrature points of one cell to + * quadrature points of another cell as well as between processors + * upon mesh refinement and repartitioning. + * The transfer of quadrature point data between parent and child cells requires + * some kind of projection and/or interpolation. + * One possible implementation is via the L2 projection and prolongation matrices + * as implemented in ContinuousQuadratureDataTransfer class. + * + * @author Denis Davydov, Jean-Paul Pelteret, 2016 + */ +class TransferableQuadraturePointData +{ +public: + /** + * Default constructor. + */ + TransferableQuadraturePointData(); + + /** + * Default virtual destructor. + */ + virtual ~TransferableQuadraturePointData(); + + /** + * Return the total number of values which will be + * packed/unpacked from the user's DataType class. Consequently it is also + * the size of the vectors in pack_values() and unpack_values() . + */ + virtual unsigned int number_of_values() const = 0; + + /** + * A virtual function that have to be implemented in derived classes to + * pack all data stored in the derived class into a vector @p values . + * This vector will contain all scalar and/or Tensorial data local to each + * quadrature point. + * + * @note The function will be called with @p values of size number_of_values(). + * The implementation may still have an assert to check that it is indeed the case. + */ + virtual void pack_values(std::vector &values) const = 0; + + /** + * The opposite of the above, namely + * unpack a vector @p values into the data stored in this class. + * + * @note The function will be called with @p values of size number_of_values(). + * The implementation may still have an assert to check that it is indeed the case. + */ + virtual void unpack_values(const std::vector &values) = 0; +}; + + +#ifdef DEAL_II_WITH_P4EST +namespace parallel +{ + namespace distributed + { + /** + * A class for the transfer of continuous data stored at quadrature points when + * performing h-adaptive refinement of parallel::distributed::Triangulation . + * + *

    Implementation details

    + * + * This class implements the transfer of the quadrature point data between + * cells in case of adaptive refinement using L2 projection. That also includes + * automatic shipping of information between different processors. + * + * To that end, the constructor of the class is provided with three main objects: + * scalar FiniteElement @projection_fe, @p mass_quadrature and @p data_quadrature + * Quadrature rules. + * First, the data located at @p data_quadrature of each cell is L2-projected + * to the continuous space defined by a single FiniteElement @projection_fe . + * This is achieved using FETools::compute_projection_from_quadrature_points_matrix(). + * In doing so the mass matrix of this element is required, which will be calculated + * with the @p mass_quadrature rule . Should the cell now belong to another processor, + * the data is then sent to this processor. The class makes use of a feature of p4est + * (and parallel::distributed::Triangulation) that allows one to attach information + * to cells during mesh refinement and rebalancing. + * On receiving information on the target cell, the data is projected back to the + * quadrature points using the matrix calculated by + * FETools::compute_interpolation_to_quadrature_points_matrix() . + * In the case that local refinement is performed, this class first + * project local DoF values of the parent element to each child. + * + * + * This class is templated by @p DataType type, however the user's @p DataType class + * has to be derived from the TransferableQuadraturePointData class. In practice + * that amounts to implementing the following three functions shown below for a + * quadrature point data with 2 scalars: + * @code + * class MyQData : public TransferableQuadraturePointData + * { + * public: + * double elasticity_parameter_lambda; + * double elasticity_parameter_mu; + * unsigned int number_of_values() const + * { + * return 2; + * } + * + * // a function to pack scalars into a vector + * void pack_values(std::vector &values) const + * { + * Assert (scalars.size()==2, ExcInternalError()); + * scalars[0] = elasticity_parameter_lambda; + * scalars[1] = elasticity_parameter_mu; + * } + * + * void unpack_values(const std::vector &values) + * { + * Assert (scalars.size() ==2, ExcInternalError()); + * elasticity_parameter_lambda = scalars[0]; + * elasticity_parameter_mu = scalars[1]; + * } + * }; + * @endcode + * Note that the order of packing and unpacking has to be the same. + * + * This class can then be use with CellDataStorage in the following way: + * @code + * CellDataStorage::cell_iterator,MyQData> data_storage; + * parallel::distributed::ContinuousQuadratureDataTransfer data_transfer(FE_Q(2),QGauss(3),QGauss(4)); + * //...populate data for all active cells in data_storage + * //...mark cells for refinement... + * data_transfer.prepare_for_coarsening_and_refinement(triangulation,data_storage); + * triangulation.execute_coarsening_and_refinement(); + * //...initialize quadrature point data on new cells by calling CellDataStorage::reinit() + * data_transfer.interpolate(); + * @endcode + * This approach can be extended to quadrature point data with Tensors of + * arbitrary order, although with a little bit more work in packing and + * unpacking of data inside MyQData class. + * + * @note Currently coarsening is not supported. + * + * @note The functionality provided by this class can alternatively be achieved + * using parallel::distributed::SolutionTransfer. However, that would require + * the following steps: (i) create an auxiliary DoFHandler with a (discontinuous Galerkin) + * FiniteElement which has enough components to represent all data stored at the + * quadrature points; (ii) project the data to the FiniteElement space + * and thereby store results in global vectors; (iii) use parallel::distributed::SolutionTransfer + * to project FE vectors to the new mesh; and (iv) finally project the + * data back to the quadrature points on the new mesh via FEValues class. + * The ContinuousQuadratureDataTransfer class aims at simplifying the whole + * process by only requiring that the quadrature point data class is derived + * from the TransferableQuadraturePointData. Everything else will be done automatically. + * + * @note This class is not well suited to situations where the values stored + * at quadrature points represent samples from a discontinuous field. An + * example for such a situation would be where the data stored at the + * quadrature points represents the elastic or plastic state of a material, + * i.e., a property that varies discontinuously within a solid. In such cases, + * trying to transfer data from the quadrature points to a finite element + * field that is continuous (at least within one cell) will likely yield over + * and undershoots that, once evaluated at a different set of quadrature points + * (on child or parent cells) results in values that will not make much sense. + * + * @author Denis Davydov, Jean-Paul Pelteret, 2016 + */ + template + class ContinuousQuadratureDataTransfer + { + public: + static_assert(std::is_base_of::value, + "User's DataType class should be derived from QPData"); + + /** + * A typedef for a cell. + */ + typedef typename parallel::distributed::Triangulation::cell_iterator CellIteratorType; + + /** + * Constructor which takes the FiniteElement @p projection_fe , the quadrature + * rule @p mass_quadrature used to integrate its local mass matrix and + * finally the quadrature rule @p data_quadrature which is used to store @p DataType. + * + * @pre @p projection_fe has to be scalar-valued. + * + * @note Since this class does projection on cell-by-cell basis, + * @p projection_fe is only required to be continuous within the cell. + */ + ContinuousQuadratureDataTransfer(const FiniteElement &projection_fe, + const Quadrature &mass_quadrature, + const Quadrature &data_quadrature); + + /** + * Destructor. + */ + ~ContinuousQuadratureDataTransfer(); + + /** + * Prepare for coarsening and refinement of a triangulation @p tria . + * @p data_storage represents the cell data which should be transferred and + * it should be initialized for each locally owned active cell. + * + * @note Although CellDataStorage class allows storing on different cells + * different objects derived from the base class, here we assume that + * @p data_storage contains objects of the same type, more specifically + * they pack/unpack the same data. + */ + void prepare_for_coarsening_and_refinement (parallel::distributed::Triangulation &tria, + CellDataStorage &data_storage); + + /** + * Interpolate the data previously stored in this object before the mesh + * was refined or coarsened onto the quadrature points of the currently active + * set of cells. + * + * @note Before calling this function the user is expected to populate the + * data stored in the @p data_storage object provided to prepare_for_coarsening_and_refinement() + * at new cells using CellDataStorage::initialize(). If that is not the case, + * an exception will be thrown in debug mode. + */ + void interpolate (); + + private: + /** + * A callback function used to pack the data on the current mesh into + * objects that can later be retrieved after refinement, coarsening and + * repartitioning. + */ + void pack_function(const typename parallel::distributed::Triangulation::cell_iterator &cell, + const typename parallel::distributed::Triangulation::CellStatus status, + void *data); + + /** + * A callback function used to unpack the data on the current mesh that + * has been packed up previously on the mesh before refinement, + * coarsening and repartitioning. + */ + void unpack_function (const typename parallel::distributed::Triangulation::cell_iterator &cell, + const typename parallel::distributed::Triangulation::CellStatus status, + const void *data); + + /** + * FiniteElement used to project data from and to quadrature points. + */ + const std::unique_ptr > projection_fe; + + /** + * The size of the data that will be sent, which depends on the DataType class. + */ + std::size_t data_size_in_bytes; + + /** + * Number of quadrature points at which DataType is stored. + */ + const unsigned int n_q_points; + + /** + * Projection matrix from the quadrature points to local DoFs for a single + * scalar. + */ + FullMatrix project_to_fe_matrix; + + /** + * Projection matrix from the local DoFs to quadrature points for a single + * scalar. + */ + FullMatrix project_to_qp_matrix; + + /** + * Auxiliary matrix which represents projection of each internal value stored + * at the quadrature point (second index) to the local DoFs of the @p projection_fe + * (first index). + */ + FullMatrix matrix_dofs; + + /** + * Projection of @p matrix_dofs to each child cell in case of adaptive refinement. + */ + FullMatrix matrix_dofs_child; + + /** + * Auxiliary matrix which represents data (second index) stored at each + * quadrature point (first index). + */ + FullMatrix matrix_quadrature; + + /** + * The offset that the parallel::distributed::Triangulation has assigned + * to this object starting at which we are allowed to read. + */ + unsigned int offset; + + /** + * A pointer to the CellDataStorage class whose data will be transfered. + */ + CellDataStorage *data_storage; + + /** + * A pointer to the distributed triangulation to which cell data is attached. + */ + parallel::distributed::Triangulation *triangulation; + }; + + } // distributed + +} // parallel + +#endif + +/*@}*/ + +#ifndef DOXYGEN + +// ------------------- inline and template functions ---------------- + +//-------------------------------------------------------------------- +// TransferableQuadraturePointData +//-------------------------------------------------------------------- +inline TransferableQuadraturePointData::TransferableQuadraturePointData() +{}; + +inline TransferableQuadraturePointData::~TransferableQuadraturePointData() +{}; + +//-------------------------------------------------------------------- +// CellDataStorage +//-------------------------------------------------------------------- + +template +CellDataStorage:: +CellDataStorage() +{}; + + + +template +CellDataStorage:: +~CellDataStorage() +{}; + + + +template +template +void CellDataStorage::initialize(const CellIteratorType &cell, + const unsigned int n_q_points) +{ + static_assert(std::is_base_of::value, + "User's T class should be derived from user's DataType class"); + static_assert(std::is_default_constructible::value, + "Class T requires default-constructible elements"); + + if (map.find(cell) == map.end()) + { + map[cell] = std::vector >(n_q_points); + // we need to initialise one-by-one as the std::vector<>(q, T()) + // will end with a single same T object stored in each element of the vector: + typename std::map >::iterator it = map.find(cell); + for (unsigned int q=0; q < n_q_points; q++) + it->second[q] = std::make_shared(); + } +} + + + +template +template +void CellDataStorage::initialize(const CellIteratorType &cell_start, + const CellIteratorType &cell_end, + const unsigned int number) +{ + for (CellIteratorType it = cell_start; it != cell_end; it++) + if (it->is_locally_owned()) + initialize(it,number); +} + + + +template +bool CellDataStorage::erase(const CellIteratorType &cell) +{ + const typename std::map >::const_iterator it = map.find(cell); + for (unsigned int i = 0; i < it->second.size(); i++) + { + Assert(it->second[i].unique(), + ExcMessage("Can not erase the cell data multiple objects reference its data.")); + } + + return (map.erase(cell)==1); +} + + + + +template +void CellDataStorage::clear() +{ + // Do not call + // map.clear(); + // as we want to be sure no one uses the stored objects. Loop manually: + typename std::map >::iterator it = map.begin(); + while (it != map.end()) + { + // loop over all objects and see if noone is using them + for (unsigned int i = 0; i < it->second.size(); i++) + { + Assert(it->second[i].unique(), + ExcMessage("Can not erase the cell data, multiple objects reference it.")); + } + it = map.erase(it); + } +} + + + + +template +template +CellDataStorage::CellVector +CellDataStorage::get_data(const CellIteratorType &cell) +{ + static_assert(std::is_base_of::value, + "User's T class should be derived from user's DataType class"); + + typename std::map >::iterator it = map.find(cell); + Assert(it != map.end(), ExcMessage("Could not find data for the cell")); + + // It would be nice to have a specialized version of this function for T==DataType. + // However explicit (i.e full) specialization of a member template is only allowed + // when the enclosing class is also explicitly (i.e fully) specialized. + // Thus, stick with copying of shared pointers even when the T==DataType: + CellVector res (it->second.size()); + for (unsigned int q = 0; q < res.size(); q++) + res[q] = std::dynamic_pointer_cast(it->second[q]); + return res; +} + + + + +template +template +CellDataStorage::CellVector +CellDataStorage::get_data(const CellIteratorType &cell) const +{ + static_assert(std::is_base_of::value, + "User's T class should be derived from user's DataType class"); + + const typename std::map >::const_iterator it = map.find(cell); + Assert(it != map.end(), ExcMessage("Could not find QP data for the cell")); + + // Cast base class to the desired class. This has to be done irrespectively of + // T==DataType as we need to return shapred_ptr to make sure the user + // does not modify the content of QP objects + CellVector res (it->second.size()); + for (unsigned int q = 0; q < res.size(); q++) + res[q] = std::dynamic_pointer_cast(it->second[q]); + return res; +} + +//-------------------------------------------------------------------- +// ContinuousQuadratureDataTransfer +//-------------------------------------------------------------------- + + +/* + * Pack cell data of type @p DataType stored using @p data_storage in @p cell + * at each quadrature point to @p matrix_data. Here @p matrix_data is a matrix + * whose first index corresponds to different quadrature points on the cell + * whereas the second index represents different values stored at each + * quadrature point in the DataType class. + */ +template +void pack_cell_data +(const CellIteratorType &cell, + const CellDataStorage *data_storage, + FullMatrix &matrix_data) +{ + static_assert(std::is_base_of::value, + "User's DataType class should be derived from QPData"); + + const std::vector > qpd = data_storage->get_data(cell); + + const unsigned int m = matrix_data.m(); + const unsigned int n = matrix_data.n(); + + std::vector single_qp_data(n); + Assert (qpd.size() == m, + ExcDimensionMismatch(qpd.size(),m)); + for (unsigned int q = 0; q < qpd.size(); q++) + { + qpd[q]->pack_values(single_qp_data); + Assert (single_qp_data.size() == n, + ExcDimensionMismatch(single_qp_data.size(),n)); + + for (unsigned int i = 0; i < n; i++) + matrix_data(q,i) = single_qp_data[i]; + } +} + + + + +/* + * the opposite of the pack function above. + */ +template +void unpack_to_cell_data +(const CellIteratorType &cell, + const FullMatrix &values_at_qp, + CellDataStorage *data_storage) +{ + static_assert(std::is_base_of::value, + "User's DataType class should be derived from QPData"); + + std::vector > qpd = data_storage->get_data(cell); + + const unsigned int m = values_at_qp.m(); + const unsigned int n = values_at_qp.n(); + + std::vector single_qp_data(n); + Assert(qpd.size() == m, + ExcDimensionMismatch(qpd.size(),m)); + + for (unsigned int q = 0; q < qpd.size(); q++) + { + for (unsigned int i = 0; i < n; i++) + single_qp_data[i] = values_at_qp(q,i); + qpd[q]->unpack_values(single_qp_data); + } +} + + +#ifdef DEAL_II_WITH_P4EST + +namespace parallel +{ + namespace distributed + { + template + ContinuousQuadratureDataTransfer::ContinuousQuadratureDataTransfer + (const FiniteElement &projection_fe_, + const Quadrature &lhs_quadrature, + const Quadrature &rhs_quadrature) + : + projection_fe(std::unique_ptr >(projection_fe_.clone())), + data_size_in_bytes(0), + n_q_points(rhs_quadrature.size()), + project_to_fe_matrix(projection_fe->dofs_per_cell,n_q_points), + project_to_qp_matrix(n_q_points,projection_fe->dofs_per_cell), + offset(0), + data_storage(nullptr), + triangulation(nullptr) + { + Assert(projection_fe->n_components() == 1, + ExcMessage("ContinuousQuadratureDataTransfer requires scalar FiniteElement")); + + FETools::compute_projection_from_quadrature_points_matrix + (*projection_fe.get(), + lhs_quadrature, + rhs_quadrature, + project_to_fe_matrix); + + FETools::compute_interpolation_to_quadrature_points_matrix + (*projection_fe.get(), + rhs_quadrature, + project_to_qp_matrix); + } + + + + template + ContinuousQuadratureDataTransfer::~ContinuousQuadratureDataTransfer() + {} + + + + template + void ContinuousQuadratureDataTransfer::prepare_for_coarsening_and_refinement + (parallel::distributed::Triangulation &tr_, + CellDataStorage &data_storage_) + { + Assert (data_storage == nullptr, + ExcMessage("This function can be called only once")); + triangulation = &tr_; + data_storage = &data_storage_; + // get the number from the first active cell + unsigned int number_of_values = 0; + // if triangulation has some active cells locally owned cells on this processor we can expect + // data to be initialized. Do that to get the number: + for (typename parallel::distributed::Triangulation::cell_iterator it = triangulation->begin_active(); + it != triangulation->end(); it++) + if (it->is_locally_owned()) + { + std::vector > qpd = data_storage->get_data(it); + number_of_values = qpd[0]->number_of_values(); + break; + } + // some processors may have no data stored, thus get the maximum among all processors: + number_of_values = Utilities::MPI::max(number_of_values, triangulation->get_communicator ()); + Assert (number_of_values > 0, + ExcInternalError()); + const unsigned int dofs_per_cell = projection_fe->dofs_per_cell; + matrix_dofs.reinit(dofs_per_cell,number_of_values); + matrix_dofs_child.reinit(dofs_per_cell,number_of_values); + matrix_quadrature.reinit(n_q_points,number_of_values); + data_size_in_bytes = sizeof(double) * dofs_per_cell * number_of_values; + + offset = triangulation->register_data_attach(data_size_in_bytes, + std_cxx11::bind(&ContinuousQuadratureDataTransfer::pack_function, + this, + std_cxx11::_1, + std_cxx11::_2, + std_cxx11::_3)); + } + + + + template + void ContinuousQuadratureDataTransfer::interpolate () + { + triangulation->notify_ready_to_unpack(offset, + std_cxx11::bind(&ContinuousQuadratureDataTransfer::unpack_function, + this, + std_cxx11::_1, + std_cxx11::_2, + std_cxx11::_3)); + + // invalidate the pointers + data_storage = nullptr; + triangulation = nullptr; + } + + + + template + void ContinuousQuadratureDataTransfer::pack_function + (const typename parallel::distributed::Triangulation::cell_iterator &cell, + const typename parallel::distributed::Triangulation::CellStatus status, + void *data) + { + double *data_store = reinterpret_cast(data); + + pack_cell_data(cell,data_storage,matrix_quadrature); + + // project to FE + project_to_fe_matrix.mmult(matrix_dofs, matrix_quadrature); + + std::memcpy(data_store, &matrix_dofs(0,0), data_size_in_bytes); + } + + + + template + void ContinuousQuadratureDataTransfer:: + unpack_function (const typename parallel::distributed::Triangulation::cell_iterator &cell, + const typename parallel::distributed::Triangulation::CellStatus status, + const void *data) + { + Assert ((status!=parallel::distributed::Triangulation::CELL_COARSEN), + ExcNotImplemented()); + + const double *data_store = reinterpret_cast(data); + + std::memcpy(&matrix_dofs(0,0), data_store, data_size_in_bytes); + + if (cell->has_children()) + { + // we need to first use prolongation matrix to get dofvalues on child + // cells based on dofvalues stored in the parent's data_store + for (unsigned int child=0; childn_children(); ++child) + if (cell->child(child)->is_locally_owned()) + { + projection_fe->get_prolongation_matrix(child, cell->refinement_case()) + .mmult (matrix_dofs_child, matrix_dofs); + + // now we do the usual business of evaluating FE on quadrature points: + project_to_qp_matrix.mmult(matrix_quadrature,matrix_dofs_child); + + // finally, put back into the map: + unpack_to_cell_data(cell->child(child),matrix_quadrature,data_storage); + } + } + else + { + // if there are no children, evaluate FE field at + // rhs_quadrature points. + project_to_qp_matrix.mmult(matrix_quadrature,matrix_dofs); + + // finally, put back into the map: + unpack_to_cell_data(cell,matrix_quadrature,data_storage); + } + } + + } // distributed + +} // parallel + +#endif // DEAL_II_WITH_P4EST + +#endif // DOXYGEN +DEAL_II_NAMESPACE_CLOSE + +#endif // CXX11 + +#endif diff --git a/tests/base/cell_data_storage_01.cc b/tests/base/cell_data_storage_01.cc new file mode 100644 index 0000000000..3e4ad199ca --- /dev/null +++ b/tests/base/cell_data_storage_01.cc @@ -0,0 +1,171 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2016 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + + + +// Check CellDataStorage initialize(cell,number), initialize(start,end,number) and get_data() functions. + + +#include "../tests.h" + +#ifdef DEAL_II_WITH_CXX11 + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include + +#include + +using namespace dealii; + +template +class MyFunction : public Function +{ +public: + MyFunction() + : + Function(1) + {} + + double value(const Point &p, + const unsigned int comp=0) const + { + const double x = p[0]; + const double y = p[1]; + // some function we know we can project with FE_Q(2) + return 0.5*x*x+2.1*y*y+2; + } +}; + +class MyQData +{ +public: + MyQData() {}; + virtual ~MyQData() {}; + + double value; +}; + +const double eps = 1e-10; +DeclException3 (ExcWrongValue, + double, double, double, + << arg1 <<" != "<< arg2 << " with delta = " << arg3 ); + + +/** + * Loop over quadrature points and check that value is the same as given by the function. + */ +template +void check_qph(Triangulation &tr, + const CellDataStorage::cell_iterator,DATA> &manager, + const Quadrature &rhs_quadrature, + const MyFunction &func) +{ + DoFHandler dof_handler(tr); + FE_Q dummy_fe(1); + FEValues fe_values(dummy_fe, rhs_quadrature, update_quadrature_points); + dof_handler.distribute_dofs( dummy_fe); + typename Triangulation::active_cell_iterator cell; + for (cell = tr.begin_active(); + cell != tr.end(); + ++cell) + if (cell->is_locally_owned()) + { + typename DoFHandler::active_cell_iterator dof_cell(*cell, &dof_handler); + fe_values.reinit(dof_cell); + const std::vector > &q_points = fe_values.get_quadrature_points(); + const std::vector > qpd = manager.get_data(cell); + for (unsigned int q=0; q < q_points.size(); q++) + { + const double value = func.value(q_points[q]); + const double value2= qpd[q]->value; + AssertThrow(std::fabs(value-value2) < eps, + ExcWrongValue(value,value2,value-value2)); + } + } + dof_handler.clear(); +} + +template +void test() +{ + const MyFunction my_func; + Triangulation tr; + + GridGenerator::subdivided_hyper_cube(tr, 2); + tr.refine_global(1); + typename Triangulation::active_cell_iterator cell; + + // pppulate quadrature point data + QGauss rhs(4); + CellDataStorage::cell_iterator,MyQData> data_storage; + { + DoFHandler dof_handler(tr); + FE_Q dummy_fe(1); + FEValues fe_values(dummy_fe, rhs, update_quadrature_points); + dof_handler.distribute_dofs( dummy_fe); + for (cell = tr.begin_active(); + cell != tr.end(); + ++cell) + if (cell->is_locally_owned()) + { + typename DoFHandler::active_cell_iterator dof_cell(*cell, &dof_handler); + fe_values.reinit(dof_cell); + const std::vector > &q_points = fe_values.get_quadrature_points(); + data_storage.initialize(cell,rhs.size()); + std::vector > qpd = data_storage.get_data(cell); + for (unsigned int q=0; q < rhs.size(); q++) + qpd[q]->value = my_func.value(q_points[q]); + } + dof_handler.clear(); + } + + check_qph(tr,data_storage,rhs,my_func); + + data_storage.initialize(tr.begin_active(),tr.end(),rhs.size()); + + check_qph(tr,data_storage,rhs,my_func); + + deallog << "Ok" << std::endl; + +} + + +int main(int argc, char *argv[]) +{ + Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1); + mpi_initlog(); + + test<2>(); +} + +#else +int main(int argc, char *argv[]) +{ + Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1); + mpi_initlog(); + deallog << "Ok" << std::endl; +} +#endif diff --git a/tests/base/cell_data_storage_01.output b/tests/base/cell_data_storage_01.output new file mode 100644 index 0000000000..6ea1f897a5 --- /dev/null +++ b/tests/base/cell_data_storage_01.output @@ -0,0 +1,2 @@ + +DEAL::Ok diff --git a/tests/base/cell_data_storage_02.cc b/tests/base/cell_data_storage_02.cc new file mode 100644 index 0000000000..9962a78c55 --- /dev/null +++ b/tests/base/cell_data_storage_02.cc @@ -0,0 +1,152 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2016 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + + + +// Check that CellDataStorage erase() and clear() + + +#include "../tests.h" + +#ifdef DEAL_II_WITH_CXX11 + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include + +#include + +using namespace dealii; + +template +class MyFunction : public Function +{ +public: + MyFunction() + : + Function(1) + {} + + double value(const Point &p, + const unsigned int comp=0) const + { + const double x = p[0]; + const double y = p[1]; + // some function we know we can project with FE_Q(2) + return 0.5*x*x+2.1*y*y+2; + } +}; + +const double default_value = 0.; + +class MyQData +{ +public: + MyQData(): value(default_value) {}; + virtual ~MyQData() {}; + + double value; +}; + +DeclException3 (ExcWrongValue, + double, double, double, + << arg1 <<" != "<< arg2 << " with delta = " << arg3 ); + + +template +void test() +{ + const MyFunction my_func; + + Triangulation tr; + + GridGenerator::subdivided_hyper_cube(tr, 2); + tr.refine_global(1); + typename Triangulation::active_cell_iterator cell; + + // pppulate quadrature point data + QGauss rhs(4); + CellDataStorage::cell_iterator,MyQData> data_storage; + { + DoFHandler dof_handler(tr); + FE_Q dummy_fe(1); + FEValues fe_values(dummy_fe, rhs, update_quadrature_points); + dof_handler.distribute_dofs( dummy_fe); + for (cell = tr.begin_active(); + cell != tr.end(); + ++cell) + if (cell->is_locally_owned()) + { + typename DoFHandler::active_cell_iterator dof_cell(*cell, &dof_handler); + fe_values.reinit(dof_cell); + const std::vector > &q_points = fe_values.get_quadrature_points(); + data_storage.initialize(cell,rhs.size()); + { + std::vector > qpd = data_storage.get_data(cell); + for (unsigned int q=0; q < rhs.size(); q++) + qpd[q]->value = my_func.value(q_points[q]); + } + + // do erase + const bool erased = data_storage.erase(cell); + Assert (erased, + ExcInternalError()); + // initialize with default constructor + data_storage.initialize(cell,rhs.size()); + // check that values are now zero (see default constructor) + { + std::vector > qpd = data_storage.get_data(cell); + for (unsigned int q=0; q < rhs.size(); q++) + AssertThrow(qpd[q]->value == default_value, + ExcWrongValue(qpd[q]->value,default_value,(qpd[q]->value-default_value))); + } + + } + dof_handler.clear(); + } + + // call clear + data_storage.clear(); + + deallog << "Ok" << std::endl; + +} + + +int main(int argc, char *argv[]) +{ + Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1); + mpi_initlog(); + + test<2>(); +} + +#else +int main(int argc, char *argv[]) +{ + Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1); + mpi_initlog(); + deallog << "Ok" << std::endl; +} +#endif diff --git a/tests/base/cell_data_storage_02.output b/tests/base/cell_data_storage_02.output new file mode 100644 index 0000000000..6ea1f897a5 --- /dev/null +++ b/tests/base/cell_data_storage_02.output @@ -0,0 +1,2 @@ + +DEAL::Ok diff --git a/tests/base/quadrature_point_data.cc b/tests/base/quadrature_point_data.cc new file mode 100644 index 0000000000..bb9cdea24a --- /dev/null +++ b/tests/base/quadrature_point_data.cc @@ -0,0 +1,219 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2016 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + + + +// Check that QuadraturePointManager class. To that end first evaluate some +// quadratic function at quadrature points. Then refine cells and project using +// FE_Q(2). Finally check that the values at quadrature points are still consistent +// with the original function. + + +#include "../tests.h" + +#ifdef DEAL_II_WITH_CXX11 + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include + +#include + +using namespace dealii; + +template +class MyFunction : public Function +{ +public: + MyFunction() + : + Function(1) + {} + + double value(const Point &p, + const unsigned int comp=0) const + { + const double x = p[0]; + const double y = p[1]; + // some function we know we can project with FE_Q(2) + return 0.5*x*x+2.1*y*y+2; + } +}; + +class MyQData : public TransferableQuadraturePointData +{ +public: + MyQData() {}; + virtual ~MyQData() {}; + + unsigned int number_of_values() const + { + return 1; + } + + double value; + + virtual void pack_values(std::vector &scalars) const + { + Assert (scalars.size()==1, ExcInternalError()); + scalars[0] = value; + } + + virtual void unpack_values(const std::vector &scalars) + { + Assert (scalars.size() ==1, ExcInternalError()); + value = scalars[0]; + } +}; + +const double eps = 1e-10; +DeclException3 (ExcWrongValue, + double, double, double, + << arg1 <<" != "<< arg2 << " with delta = " << arg3 ); + + +/** + * Loop over quadrature points and check that value is the same as given by the function. + */ +template +void check_qph(parallel::distributed::Triangulation &tr, + const CellDataStorage::cell_iterator,DATA> &manager, + const Quadrature &rhs_quadrature, + const MyFunction &func) +{ + DoFHandler dof_handler(tr); + FE_Q dummy_fe(1); + FEValues fe_values(dummy_fe, rhs_quadrature, update_quadrature_points); + dof_handler.distribute_dofs( dummy_fe); + typename Triangulation::active_cell_iterator cell; + for (cell = tr.begin_active(); + cell != tr.end(); + ++cell) + if (cell->is_locally_owned()) + { + typename DoFHandler::active_cell_iterator dof_cell(*cell, &dof_handler); + fe_values.reinit(dof_cell); + const std::vector > &q_points = fe_values.get_quadrature_points(); + const std::vector > qpd = manager.get_data(cell); + for (unsigned int q=0; q < q_points.size(); q++) + { + const double value = func.value(q_points[q]); + const double value2= qpd[q]->value; + AssertThrow(std::fabs(value-value2) < eps, + ExcWrongValue(value,value2,value-value2)); + } + } + dof_handler.clear(); +} + +template +void test() +{ + unsigned int myid = Utilities::MPI::this_mpi_process (MPI_COMM_WORLD); + unsigned int numprocs = Utilities::MPI::n_mpi_processes (MPI_COMM_WORLD); + + const MyFunction my_func; + + parallel::distributed::Triangulation tr(MPI_COMM_WORLD); + + GridGenerator::subdivided_hyper_cube(tr, 2); + tr.refine_global(1); + typename Triangulation::active_cell_iterator cell; + + // pppulate quadrature point data + QGauss rhs(4); + CellDataStorage::cell_iterator,MyQData> data_storage; + parallel::distributed::ContinuousQuadratureDataTransfer data_transfer(FE_Q(2),QGauss(3),rhs); + { + DoFHandler dof_handler(tr); + FE_Q dummy_fe(1); + FEValues fe_values(dummy_fe, rhs, update_quadrature_points); + dof_handler.distribute_dofs( dummy_fe); + for (cell = tr.begin_active(); + cell != tr.end(); + ++cell) + if (cell->is_locally_owned()) + { + typename DoFHandler::active_cell_iterator dof_cell(*cell, &dof_handler); + fe_values.reinit(dof_cell); + const std::vector > &q_points = fe_values.get_quadrature_points(); + data_storage.initialize(cell,rhs.size()); + std::vector > qpd = data_storage.get_data(cell); + for (unsigned int q=0; q < rhs.size(); q++) + qpd[q]->value = my_func.value(q_points[q]); + } + dof_handler.clear(); + } + + check_qph(tr,data_storage,rhs,my_func); + + + // mark some for refinement + for (cell = tr.begin_active(); + cell != tr.end(); + ++cell) + if (cell->center()[0] < 0.5) + cell->set_refine_flag(); + + data_transfer.prepare_for_coarsening_and_refinement(tr,data_storage); + + tr.execute_coarsening_and_refinement(); + + // create qhp data + for (cell = tr.begin_active(); + cell != tr.end(); + ++cell) + if (cell->is_locally_owned()) + { + data_storage.initialize(cell,rhs.size()); + } + + data_transfer.interpolate(); + + // check that projected data still consistent with what we expect, i.e. + // exact reproduction of the function + check_qph(tr,data_storage,rhs,my_func); + + deallog << "Ok" << std::endl; + +} + + +int main(int argc, char *argv[]) +{ + Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1); + mpi_initlog(); + + test<2>(); +} + +#else +int main(int argc, char *argv[]) +{ + Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1); + mpi_initlog(); + deallog << "Ok" << std::endl; +} +#endif diff --git a/tests/base/quadrature_point_data.with_p4est=on.mpirun=3.output b/tests/base/quadrature_point_data.with_p4est=on.mpirun=3.output new file mode 100644 index 0000000000..6ea1f897a5 --- /dev/null +++ b/tests/base/quadrature_point_data.with_p4est=on.mpirun=3.output @@ -0,0 +1,2 @@ + +DEAL::Ok diff --git a/tests/base/quadrature_point_data_02.cc b/tests/base/quadrature_point_data_02.cc new file mode 100644 index 0000000000..2dbb269c7a --- /dev/null +++ b/tests/base/quadrature_point_data_02.cc @@ -0,0 +1,233 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2016 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + + + +// Check that QuadraturePointManager class. To that end first evaluate some +// quadratic function at quadrature points. Then refine cells and project using +// FE_Q(2). Finally check that the values at quadrature points are still consistent +// with the original function. +// Same as quadrature_point_data.cc, but for DATA class with 2 scalars. + + +#include "../tests.h" + +#ifdef DEAL_II_WITH_CXX11 +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include + +#include + +using namespace dealii; + +template +class MyFunction : public Function +{ +public: + MyFunction() + : + Function(1) + {} + + double value(const Point &p, + const unsigned int comp) const + { + const double x = p[0]; + const double y = p[1]; + // some function we know we can project with FE_Q(2) + if (comp ==0) + return 0.5*x*x+2.1*y*y+2; + else + return 0.1*x + 22.5*y*y; + } +}; + +struct MyQData : public TransferableQuadraturePointData +{ + double value1; + double value2; + + unsigned int number_of_values() const + { + return 2; + } + + + virtual void pack_values(std::vector &scalars) const + { + Assert (scalars.size()==2, ExcInternalError()); + scalars[0] = value1; + scalars[1] = value2; + } + + virtual void unpack_values(const std::vector &scalars) + { + Assert (scalars.size() ==2, ExcInternalError()); + value1 = scalars[0]; + value2 = scalars[1]; + } +}; + +const double eps = 1e-10; +DeclException3 (ExcWrongValue, + double, double, double, + << arg1 <<" != "<< arg2 << " with delta = " << arg3 ); + + +/** + * Loop over quadrature points and check that value is the same as given by the function. + */ +template +void check_qph(parallel::distributed::Triangulation &tr, + const CellDataStorage::cell_iterator,DATA> &manager, + const Quadrature &rhs_quadrature, + const MyFunction &func) +{ + DoFHandler dof_handler(tr); + FE_Q dummy_fe(1); + FEValues fe_values(dummy_fe, rhs_quadrature, update_quadrature_points); + dof_handler.distribute_dofs( dummy_fe); + typename Triangulation::active_cell_iterator cell; + for (cell = tr.begin_active(); + cell != tr.end(); + ++cell) + if (cell->is_locally_owned()) + { + typename DoFHandler::active_cell_iterator dof_cell(*cell, &dof_handler); + fe_values.reinit(dof_cell); + const std::vector > &q_points = fe_values.get_quadrature_points(); + const std::vector > qpd = manager.get_data(cell); + for (unsigned int q=0; q < q_points.size(); q++) + { + const double f_1 = func.value(q_points[q],0); + const double f_2 = func.value(q_points[q],1); + const double d_1= qpd[q]->value1; + const double d_2= qpd[q]->value2; + AssertThrow(std::fabs(f_1-d_1) < eps, + ExcWrongValue(f_1,d_1,f_1-d_1)); + AssertThrow(std::fabs(f_2-d_2) < eps, + ExcWrongValue(f_2,d_2,f_2-d_2)); + } + } + dof_handler.clear(); +} + +template +void test() +{ + unsigned int myid = Utilities::MPI::this_mpi_process (MPI_COMM_WORLD); + unsigned int numprocs = Utilities::MPI::n_mpi_processes (MPI_COMM_WORLD); + + const MyFunction my_func; + + parallel::distributed::Triangulation tr(MPI_COMM_WORLD); + + GridGenerator::subdivided_hyper_cube(tr, 2); + tr.refine_global(1); + typename Triangulation::active_cell_iterator cell; + + // pppulate quadrature point data + QGauss rhs(4); + CellDataStorage::cell_iterator,MyQData> data_storage; + parallel::distributed::ContinuousQuadratureDataTransfer data_transfer(FE_Q(2),QGauss(3),rhs); + + { + DoFHandler dof_handler(tr); + FE_Q dummy_fe(1); + FEValues fe_values(dummy_fe, rhs, update_quadrature_points); + dof_handler.distribute_dofs( dummy_fe); + for (cell = tr.begin_active(); + cell != tr.end(); + ++cell) + if (cell->is_locally_owned()) + { + typename DoFHandler::active_cell_iterator dof_cell(*cell, &dof_handler); + fe_values.reinit(dof_cell); + const std::vector > &q_points = fe_values.get_quadrature_points(); + data_storage.initialize(cell,rhs.size()); + std::vector > qpd = data_storage.get_data(cell); + for (unsigned int q=0; q < rhs.size(); q++) + { + qpd[q]->value1 = my_func.value(q_points[q],0); + qpd[q]->value2 = my_func.value(q_points[q],1); + } + } + dof_handler.clear(); + } + + check_qph(tr,data_storage,rhs,my_func); + + + // mark some for refinement + for (cell = tr.begin_active(); + cell != tr.end(); + ++cell) + if (cell->center()[0] < 0.5) + cell->set_refine_flag(); + + data_transfer.prepare_for_coarsening_and_refinement(tr,data_storage); + + tr.execute_coarsening_and_refinement(); + + // create qhp data + for (cell = tr.begin_active(); + cell != tr.end(); + ++cell) + if (cell->is_locally_owned()) + { + data_storage.initialize(cell,rhs.size()); + } + + data_transfer.interpolate(); + + // check that projected data still consistent with what we expect, i.e. + // exact reproduction of the function + check_qph(tr,data_storage,rhs,my_func); + + deallog << "Ok" << std::endl; + +} + + +int main(int argc, char *argv[]) +{ + Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1); + mpi_initlog(); + + test<2>(); +} + +#else +int main(int argc, char *argv[]) +{ + Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1); + mpi_initlog(); + + deallog << "Ok" << std::endl; +} + + +#endif diff --git a/tests/base/quadrature_point_data_02.with_p4est=on.mpirun=3.output b/tests/base/quadrature_point_data_02.with_p4est=on.mpirun=3.output new file mode 100644 index 0000000000..6ea1f897a5 --- /dev/null +++ b/tests/base/quadrature_point_data_02.with_p4est=on.mpirun=3.output @@ -0,0 +1,2 @@ + +DEAL::Ok diff --git a/tests/base/quadrature_point_data_03.cc b/tests/base/quadrature_point_data_03.cc new file mode 100644 index 0000000000..4642896d61 --- /dev/null +++ b/tests/base/quadrature_point_data_03.cc @@ -0,0 +1,238 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2016 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + + + +// Same as quadrature_point_data_2.cc, but for different base and actuall data classes. + + +#include "../tests.h" + +#ifdef DEAL_II_WITH_CXX11 +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include + +#include + +using namespace dealii; + +template +class MyFunction : public Function +{ +public: + MyFunction() + : + Function(1) + {} + + double value(const Point &p, + const unsigned int comp) const + { + const double x = p[0]; + const double y = p[1]; + // some function we know we can project with FE_Q(2) + if (comp ==0) + return 0.5*x*x+2.1*y*y+2; + else + return 0.1*x + 22.5*y*y; + } +}; + +struct MyDataBase : public TransferableQuadraturePointData +{ + MyDataBase() {}; + virtual ~MyDataBase() {}; +}; + + +struct MyData : public MyDataBase +{ + double value1; + double value2; + + unsigned int number_of_values() const + { + return 2; + } + + + MyData() {}; + virtual ~MyData() {}; + + virtual void pack_values(std::vector &scalars) const + { + Assert (scalars.size()==2, ExcInternalError()); + scalars[0] = value1; + scalars[1] = value2; + } + + virtual void unpack_values(const std::vector &scalars) + { + Assert (scalars.size() ==2, ExcInternalError()); + value1 = scalars[0]; + value2 = scalars[1]; + } +}; + +const double eps = 1e-10; +DeclException3 (ExcWrongValue, + double, double, double, + << arg1 <<" != "<< arg2 << " with delta = " << arg3 ); + + +/** + * Loop over quadrature points and check that value is the same as given by the function. + */ +template +void check_qph(parallel::distributed::Triangulation &tr, + const CellDataStorage::cell_iterator,MyDataBase> &manager, + const Quadrature &rhs_quadrature, + const MyFunction &func) +{ + DoFHandler dof_handler(tr); + FE_Q dummy_fe(1); + FEValues fe_values(dummy_fe, rhs_quadrature, update_quadrature_points); + dof_handler.distribute_dofs( dummy_fe); + typename Triangulation::active_cell_iterator cell; + for (cell = tr.begin_active(); + cell != tr.end(); + ++cell) + if (cell->is_locally_owned()) + { + typename DoFHandler::active_cell_iterator dof_cell(*cell, &dof_handler); + fe_values.reinit(dof_cell); + const std::vector > &q_points = fe_values.get_quadrature_points(); + const std::vector > qpd = manager.template get_data(cell); + for (unsigned int q=0; q < q_points.size(); q++) + { + const double f_1 = func.value(q_points[q],0); + const double f_2 = func.value(q_points[q],1); + const double d_1= qpd[q]->value1; + const double d_2= qpd[q]->value2; + AssertThrow(std::fabs(f_1-d_1) < eps, + ExcWrongValue(f_1,d_1,f_1-d_1)); + AssertThrow(std::fabs(f_2-d_2) < eps, + ExcWrongValue(f_2,d_2,f_2-d_2)); + } + } + dof_handler.clear(); +} + +template +void test() +{ + unsigned int myid = Utilities::MPI::this_mpi_process (MPI_COMM_WORLD); + unsigned int numprocs = Utilities::MPI::n_mpi_processes (MPI_COMM_WORLD); + + const MyFunction my_func; + + parallel::distributed::Triangulation tr(MPI_COMM_WORLD); + + GridGenerator::subdivided_hyper_cube(tr, 2); + tr.refine_global(1); + typename Triangulation::active_cell_iterator cell; + + // pppulate quadrature point data + QGauss rhs(4); + CellDataStorage::cell_iterator,MyDataBase> data_storage; + parallel::distributed::ContinuousQuadratureDataTransfer data_transfer(FE_Q(2),QGauss(3),rhs); + { + DoFHandler dof_handler(tr); + FE_Q dummy_fe(1); + FEValues fe_values(dummy_fe, rhs, update_quadrature_points); + dof_handler.distribute_dofs( dummy_fe); + for (cell = tr.begin_active(); + cell != tr.end(); + ++cell) + if (cell->is_locally_owned()) + { + typename DoFHandler::active_cell_iterator dof_cell(*cell, &dof_handler); + fe_values.reinit(dof_cell); + const std::vector > &q_points = fe_values.get_quadrature_points(); + data_storage.template initialize(cell, rhs.size()); + std::vector > qpd = data_storage.template get_data(cell); + for (unsigned int q=0; q < rhs.size(); q++) + { + qpd[q]->value1 = my_func.value(q_points[q],0); + qpd[q]->value2 = my_func.value(q_points[q],1); + } + } + dof_handler.clear(); + } + + check_qph(tr,data_storage,rhs,my_func); + + + // mark some for refinement + for (cell = tr.begin_active(); + cell != tr.end(); + ++cell) + if (cell->center()[0] < 0.5) + cell->set_refine_flag(); + + data_transfer.prepare_for_coarsening_and_refinement(tr,data_storage); + + tr.execute_coarsening_and_refinement(); + + // create qhp data + for (cell = tr.begin_active(); + cell != tr.end(); + ++cell) + if (cell->is_locally_owned()) + { + data_storage.template initialize(cell, rhs.size()); + } + + data_transfer.interpolate(); + + // check that projected data still consistent with what we expect, i.e. + // exact reproduction of the function + check_qph(tr,data_storage,rhs,my_func); + + deallog << "Ok" << std::endl; + +} + + +int main(int argc, char *argv[]) +{ + Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1); + mpi_initlog(); + + test<2>(); +} + +#else +int main(int argc, char *argv[]) +{ + Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1); + mpi_initlog(); + + deallog << "Ok" << std::endl; +} + + +#endif diff --git a/tests/base/quadrature_point_data_03.with_p4est=on.mpirun=3.output b/tests/base/quadrature_point_data_03.with_p4est=on.mpirun=3.output new file mode 100644 index 0000000000..6ea1f897a5 --- /dev/null +++ b/tests/base/quadrature_point_data_03.with_p4est=on.mpirun=3.output @@ -0,0 +1,2 @@ + +DEAL::Ok -- 2.39.5