1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2016 - 2020 by the deal.II authors
4 //
5 // This file is part of the deal.II library.
6 //
7 // The deal.II library is free software; you can use it, redistribute
8 // it, and/or modify it under the terms of the GNU Lesser General
9 // Public License as published by the Free Software Foundation; either
10 // version 2.1 of the License, or (at your option) any later version.
11 // The full text of the license can be found in the file LICENSE.md at
12 // the top level directory of deal.II.
13 //
14 // ---------------------------------------------------------------------
15 
16 #ifndef dealii_fe_series_h
17 #define dealii_fe_series_h
18 
19 
20 
21 #include <deal.II/base/config.h>
22 
23 #include <deal.II/base/exceptions.h>
24 #include <deal.II/base/subscriptor.h>
25 #include <deal.II/base/table.h>
26 #include <deal.II/base/table_indices.h>
27 #include <deal.II/base/tensor.h>
28 
29 #include <deal.II/hp/fe_collection.h>
30 #include <deal.II/hp/q_collection.h>
31 
32 #include <deal.II/lac/full_matrix.h>
33 #include <deal.II/lac/vector.h>
34 
35 #include <deal.II/numerics/vector_tools_common.h>
36 
37 #include <memory>
38 #include <string>
39 #include <vector>
40 
41 
42 DEAL_II_NAMESPACE_OPEN
43 
44 
45 /*!@addtogroup feall */
46 /*@{*/
47 
48 
49 /**
50  * This namespace offers functions to calculate expansion series of the
51  * solution on the reference element. Coefficients of expansion are often used
52  * to estimate local smoothness of the underlying FiniteElement field to decide
53  * on h- or p-adaptive refinement strategy.
54  */
55 namespace FESeries
56 {
57   /**
58    * A class to calculate expansion of a scalar FE field into Fourier series
59    * on a reference element. The exponential form of the Fourier series is
60    * based on completeness and Hermitian orthogonality of the set of exponential
61    * functions $ \phi_{\bf k}({\bf x}) = \exp(2 \pi i\, {\bf k} \cdot {\bf x})$.
62    * For example in 1D the L2-orthogonality condition reads
63    * @f[
64    *   \int_0^1 \phi_k(x) \phi_l^\ast(x) dx=\delta_{kl}.
65    * @f]
66    * Note that $ \phi_{\bf k} = \phi_{-\bf k}^\ast $.
67    *
68    * The arbitrary scalar FE field on the reference element can be expanded in
69    * the complete orthogonal exponential basis as
70    * @f[
71    *    u({\bf x})
72    *    = \sum_{\bf k} c_{\bf k} \phi_{\bf k}({\bf x}).
73    * @f]
74    * From the orthogonality property of the basis, it follows that
75    * @f[
76    *    c_{\bf k} =
77    *    \int_{[0,1]^d} u({\bf x}) \phi_{\bf k}^\ast ({\bf x}) d{\bf x}\,.
78    * @f]
79    * It is this complex-valued expansion coefficients, that are calculated by
80    * this class. Note that $ u({\bf x}) = \sum_i u_i N_i({\bf x})$,
81    * where $ N_i({\bf x}) $ are real-valued FiniteElement shape functions.
82    * Consequently $ c_{\bf k} \equiv c_{-\bf k}^\ast $ and
83    * we only need to compute $ c_{\bf k} $ for positive indices
84    * $ \bf k $ .
85    */
86   template <int dim, int spacedim = dim>
87   class Fourier : public Subscriptor
88   {
89   public:
90     using CoefficientType = typename std::complex<double>;
91 
92     /**
93      * Constructor that initializes all required data structures.
94      *
95      * The @p n_coefficients_per_direction defines the number of coefficients in
96      * each direction, @p fe_collection is the hp::FECollection for which
97      * expansion will be used and @p q_collection is the hp::QCollection used to
98      * integrate the expansion for each FiniteElement in @p fe_collection.
99      */
100     Fourier(const std::vector<unsigned int> &      n_coefficients_per_direction,
101             const hp::FECollection<dim, spacedim> &fe_collection,
102             const hp::QCollection<dim> &           q_collection);
103 
104     /**
105      * A non-default constructor. The @p n_coefficients_per_direction defines the
106      * number of modes in each direction, @p fe_collection is the hp::FECollection
107      * for which expansion will be used and @p q_collection is the hp::QCollection
108      * used to integrate the expansion for each FiniteElement
109      * in @p fe_collection.
110      *
111      * @deprecated Use a different constructor instead.
112      */
113     DEAL_II_DEPRECATED
114     Fourier(const unsigned int                     n_coefficients_per_direction,
115             const hp::FECollection<dim, spacedim> &fe_collection,
116             const hp::QCollection<dim> &           q_collection);
117 
118     /**
119      * Calculate @p fourier_coefficients of the cell vector field given by
120      * @p local_dof_values corresponding to FiniteElement with
121      * @p cell_active_fe_index .
122      */
123     template <typename Number>
124     void
125     calculate(const dealii::Vector<Number> &local_dof_values,
126               const unsigned int            cell_active_fe_index,
127               Table<dim, CoefficientType> & fourier_coefficients);
128 
129     /**
130      * Return the number of coefficients in each coordinate direction for the
131      * finite element associated with @p index in the provided hp::FECollection.
132      */
133     unsigned int
134     get_n_coefficients_per_direction(const unsigned int index) const;
135 
136     /**
137      * Calculate all transformation matrices to transfer the finite element
138      * solution to the series expansion representation.
139      *
140      * These matrices will be generated on demand by calling calculate() and
141      * stored for recurring purposes. Usually, this operation consumes a lot of
142      * workload. With this function, all matrices will be calculated in advance.
143      * This way, we can separate their costly generation from the actual
144      * application.
145      */
146     void
147     precalculate_all_transformation_matrices();
148 
149     /**
150      * Write all transformation matrices of this object to a stream for the
151      * purpose of serialization.
152      *
153      * Since any of its transformation matrices has to be generated only once
154      * for a given scenario, it is common practice to determine them in advance
155      * calling precalculate_all_transformation_matrices() and keep them via
156      * serialization.
157      */
158     template <class Archive>
159     void
160     save_transformation_matrices(Archive &ar, const unsigned int version);
161 
162     /**
163      * Read all transformation matrices from a stream and recover them for this
164      * object.
165      */
166     template <class Archive>
167     void
168     load_transformation_matrices(Archive &ar, const unsigned int version);
169 
170     /**
171      * Test for equality of two series expansion objects.
172      */
173     bool
174     operator==(const Fourier<dim, spacedim> &fourier) const;
175 
176   private:
177     /**
178      * Number of coefficients in each direction for each finite element in the
179      * registered hp::FECollection.
180      */
181     const std::vector<unsigned int> n_coefficients_per_direction;
182 
183     /**
184      * hp::FECollection for which transformation matrices will be calculated.
185      */
186     SmartPointer<const hp::FECollection<dim, spacedim>> fe_collection;
187 
188     /**
189      * hp::QCollection used in calculation of transformation matrices.
190      */
191     const hp::QCollection<dim> q_collection;
192 
193     /**
194      * Angular frequencies $ 2 \pi {\bf k} $ .
195      */
196     Table<dim, Tensor<1, dim>> k_vectors;
197 
198     /**
199      * Transformation matrices for each FiniteElement.
200      */
201     std::vector<FullMatrix<CoefficientType>> fourier_transform_matrices;
202 
203     /**
204      * Auxiliary vector to store unrolled coefficients.
205      */
206     std::vector<CoefficientType> unrolled_coefficients;
207   };
208 
209 
210 
211   /**
212    * A class to calculate expansion of a scalar FE field into series of Legendre
213    * functions on a reference element.
214    *
215    * Legendre functions are solutions to Legendre's differential equation
216    * @f[
217    *    \frac{d}{dx}\left([1-x^2] \frac{d}{dx} P_n(x)\right) +
218    *    n[n+1] P_n(x) = 0
219    * @f]
220    * and can be expressed using Rodrigues' formula
221    * @f[
222    *    P_n(x) = \frac{1}{2^n n!} \frac{d^n}{dx^n}[x^2-1]^n.
223    * @f]
224    * These polynomials are orthogonal with respect to the $ L^2 $ inner
225    * product on the interval $ [-1;1] $
226    * @f[
227    *    \int_{-1}^1 P_m(x) P_n(x) = \frac{2}{2n + 1} \delta_{mn}
228    * @f]
229    * and are complete.
230    * A family of $ L^2 $-orthogonal polynomials on $ [0;1] $ can be
231    * constructed via
232    * @f[
233    *    \widetilde P_m = \sqrt{2} P_m(2x-1).
234    * @f]
235    *
236    *
237    * An arbitrary scalar FE field on the reference element $ [0;1] $ can be
238    * expanded in the complete orthogonal basis as
239    * @f[
240    *    u(x)
241    *    = \sum_{m} c_m \widetilde P_{m}(x).
242    * @f]
243    * From the orthogonality property of the basis, it follows that
244    * @f[
245    *    c_m = \frac{2m+1}{2}
246    *    \int_0^1 u(x) \widetilde P_m(x) dx .
247    * @f]
248    * This class calculates coefficients $ c_{\bf k} $ using
249    * $ dim $-dimensional Legendre polynomials constructed from
250    * $ \widetilde P_m(x) $ using tensor product rule.
251    */
252   template <int dim, int spacedim = dim>
253   class Legendre : public Subscriptor
254   {
255   public:
256     using CoefficientType = double;
257 
258     /**
259      * Constructor that initializes all required data structures.
260      *
261      * The @p n_coefficients_per_direction defines the number of coefficients in
262      * each direction, @p fe_collection is the hp::FECollection for which
263      * expansion will be used and @p q_collection is the hp::QCollection used to
264      * integrate the expansion for each FiniteElement in @p fe_collection.
265      */
266     Legendre(const std::vector<unsigned int> &n_coefficients_per_direction,
267              const hp::FECollection<dim, spacedim> &fe_collection,
268              const hp::QCollection<dim> &           q_collection);
269 
270     /**
271      * A non-default constructor. The @p size_in_each_direction defines the number
272      * of coefficients in each direction, @p fe_collection is the hp::FECollection
273      * for which expansion will be used and @p q_collection is the hp::QCollection
274      * used to integrate the expansion for each FiniteElement in @p fe_collection.
275      *
276      * @deprecated Use a different constructor instead.
277      */
278     DEAL_II_DEPRECATED
279     Legendre(const unsigned int n_coefficients_per_direction,
280              const hp::FECollection<dim, spacedim> &fe_collection,
281              const hp::QCollection<dim> &           q_collection);
282 
283     /**
284      * Calculate @p legendre_coefficients of the cell vector field given by
285      * @p local_dof_values corresponding to FiniteElement with
286      * @p cell_active_fe_index .
287      */
288     template <typename Number>
289     void
290     calculate(const dealii::Vector<Number> &local_dof_values,
291               const unsigned int            cell_active_fe_index,
292               Table<dim, CoefficientType> & legendre_coefficients);
293 
294     /**
295      * Return the number of coefficients in each coordinate direction for the
296      * finite element associated with @p index in the provided hp::FECollection.
297      */
298     unsigned int
299     get_n_coefficients_per_direction(const unsigned int index) const;
300 
301     /**
302      * Calculate all transformation matrices to transfer the finite element
303      * solution to the series expansion representation.
304      *
305      * These matrices will be generated on demand by calling calculate() and
306      * stored for recurring purposes. Usually, this operation consumes a lot of
307      * workload. With this function, all matrices will be calculated in advance.
308      * This way, we can separate their costly generation from the actual
309      * application.
310      */
311     void
312     precalculate_all_transformation_matrices();
313 
314     /**
315      * Write all transformation matrices of this object to a stream for the
316      * purpose of serialization.
317      *
318      * Since any of its transformation matrices has to be generated only once
319      * for a given scenario, it is common practice to determine them in advance
320      * calling precalculate_all_transformation_matrices() and keep them via
321      * serialization.
322      */
323     template <class Archive>
324     void
325     save_transformation_matrices(Archive &ar, const unsigned int version);
326 
327     /**
328      * Read all transformation matrices from a stream and recover them for this
329      * object.
330      */
331     template <class Archive>
332     void
333     load_transformation_matrices(Archive &ar, const unsigned int version);
334 
335     /**
336      * Test for equality of two series expansion objects.
337      */
338     bool
339     operator==(const Legendre<dim, spacedim> &legendre) const;
340 
341   private:
342     /**
343      * Number of coefficients in each direction for each finite element in the
344      * registered hp::FECollection.
345      */
346     const std::vector<unsigned int> n_coefficients_per_direction;
347 
348     /**
349      * hp::FECollection for which transformation matrices will be calculated.
350      */
351     SmartPointer<const hp::FECollection<dim, spacedim>> fe_collection;
352 
353     /**
354      * hp::QCollection used in calculation of transformation matrices.
355      */
356     const hp::QCollection<dim> q_collection;
357 
358     /**
359      * Transformation matrices for each FiniteElement.
360      */
361     std::vector<FullMatrix<CoefficientType>> legendre_transform_matrices;
362 
363     /**
364      * Auxiliary vector to store unrolled coefficients.
365      */
366     std::vector<CoefficientType> unrolled_coefficients;
367   };
368 
369 
370 
371   /**
372    * Calculate the @p norm of subsets of @p coefficients defined by
373    * @p predicate being constant. Return the pair of vectors of predicate values
374    * and the vector of calculated subset norms.
375    *
376    * @p predicate should return a pair of <code>bool</code> and <code>unsigned
377    * int</code>. The former is a flag whether a given TableIndices should be
378    * used in calculation, whereas the latter is the unrolled value of indices
379    * according to which the subsets of coefficients will be formed.
380    *
381    * Only those coefficients will be considered which are larger than
382    * @p smallest_abs_coefficient.
383    *
384    * @note Only the following values of @p norm_type are implemented and make
385    * sense in this case: mean, L1_norm, L2_norm, Linfty_norm. The mean norm ca
386    * only be applied to real valued coefficients.
387    */
388   template <int dim, typename CoefficientType>
389   std::pair<std::vector<unsigned int>, std::vector<double>>
390   process_coefficients(const Table<dim, CoefficientType> &coefficients,
391                        const std::function<std::pair<bool, unsigned int>(
392                          const TableIndices<dim> &)> &    predicate,
393                        const VectorTools::NormType        norm_type,
394                        const double smallest_abs_coefficient = 1e-10);
395 
396   /**
397    * Linear regression least-square fit of $y = k \, x + b$.
398    * The size of the input vectors should be equal and more than 1.
399    * The returned pair will contain $k$ (first) and $b$ (second).
400    */
401   std::pair<double, double>
402   linear_regression(const std::vector<double> &x, const std::vector<double> &y);
403 
404 } // namespace FESeries
405 
406 /*@}*/
407 
408 
409 
410 #ifndef DOXYGEN
411 
412 // -------------------  inline and template functions ----------------
413 
414 namespace internal
415 {
416   namespace FESeriesImplementation
417   {
418     template <int dim, typename CoefficientType>
419     void
fill_map_index(const Table<dim,CoefficientType> & coefficients,const TableIndices<dim> & ind,const std::function<std::pair<bool,unsigned int> (const TableIndices<dim> &)> & predicate,std::map<unsigned int,std::vector<CoefficientType>> & pred_to_values)420     fill_map_index(
421       const Table<dim, CoefficientType> &coefficients,
422       const TableIndices<dim> &          ind,
423       const std::function<
424         std::pair<bool, unsigned int>(const TableIndices<dim> &)> &predicate,
425       std::map<unsigned int, std::vector<CoefficientType>> &pred_to_values)
426     {
427       const std::pair<bool, unsigned int> pred_pair = predicate(ind);
428       // don't add a value if predicate is false
429       if (pred_pair.first == false)
430         return;
431 
432       const unsigned int     pred_value  = pred_pair.second;
433       const CoefficientType &coeff_value = coefficients(ind);
434       // If pred_value is not in the pred_to_values map, the element will be
435       // created. Otherwise a reference to the existing element is returned.
436       pred_to_values[pred_value].push_back(coeff_value);
437     }
438 
439 
440 
441     template <typename CoefficientType>
442     void
fill_map(const Table<1,CoefficientType> & coefficients,const std::function<std::pair<bool,unsigned int> (const TableIndices<1> &)> & predicate,std::map<unsigned int,std::vector<CoefficientType>> & pred_to_values)443     fill_map(
444       const Table<1, CoefficientType> &coefficients,
445       const std::function<
446         std::pair<bool, unsigned int>(const TableIndices<1> &)> &predicate,
447       std::map<unsigned int, std::vector<CoefficientType>> &     pred_to_values)
448     {
449       for (unsigned int i = 0; i < coefficients.size(0); i++)
450         {
451           const TableIndices<1> ind(i);
452           fill_map_index(coefficients, ind, predicate, pred_to_values);
453         }
454     }
455 
456 
457 
458     template <typename CoefficientType>
459     void
fill_map(const Table<2,CoefficientType> & coefficients,const std::function<std::pair<bool,unsigned int> (const TableIndices<2> &)> & predicate,std::map<unsigned int,std::vector<CoefficientType>> & pred_to_values)460     fill_map(
461       const Table<2, CoefficientType> &coefficients,
462       const std::function<
463         std::pair<bool, unsigned int>(const TableIndices<2> &)> &predicate,
464       std::map<unsigned int, std::vector<CoefficientType>> &     pred_to_values)
465     {
466       for (unsigned int i = 0; i < coefficients.size(0); i++)
467         for (unsigned int j = 0; j < coefficients.size(1); j++)
468           {
469             const TableIndices<2> ind(i, j);
470             fill_map_index(coefficients, ind, predicate, pred_to_values);
471           }
472     }
473 
474 
475 
476     template <typename CoefficientType>
477     void
fill_map(const Table<3,CoefficientType> & coefficients,const std::function<std::pair<bool,unsigned int> (const TableIndices<3> &)> & predicate,std::map<unsigned int,std::vector<CoefficientType>> & pred_to_values)478     fill_map(
479       const Table<3, CoefficientType> &coefficients,
480       const std::function<
481         std::pair<bool, unsigned int>(const TableIndices<3> &)> &predicate,
482       std::map<unsigned int, std::vector<CoefficientType>> &     pred_to_values)
483     {
484       for (unsigned int i = 0; i < coefficients.size(0); i++)
485         for (unsigned int j = 0; j < coefficients.size(1); j++)
486           for (unsigned int k = 0; k < coefficients.size(2); k++)
487             {
488               const TableIndices<3> ind(i, j, k);
489               fill_map_index(coefficients, ind, predicate, pred_to_values);
490             }
491     }
492 
493 
494 
495     template <typename Number>
496     double
complex_mean_value(const Number & value)497     complex_mean_value(const Number &value)
498     {
499       return value;
500     }
501 
502 
503 
504     template <typename Number>
505     double
complex_mean_value(const std::complex<Number> & value)506     complex_mean_value(const std::complex<Number> &value)
507     {
508       AssertThrow(false,
509                   ExcMessage(
510                     "FESeries::process_coefficients() can not be used with "
511                     "complex-valued coefficients and VectorTools::mean norm."));
512       return std::abs(value);
513     }
514   } // namespace FESeriesImplementation
515 } // namespace internal
516 
517 
518 
519 template <int dim, typename CoefficientType>
520 std::pair<std::vector<unsigned int>, std::vector<double>>
process_coefficients(const Table<dim,CoefficientType> & coefficients,const std::function<std::pair<bool,unsigned int> (const TableIndices<dim> &)> & predicate,const VectorTools::NormType norm_type,const double smallest_abs_coefficient)521 FESeries::process_coefficients(
522   const Table<dim, CoefficientType> &coefficients,
523   const std::function<std::pair<bool, unsigned int>(const TableIndices<dim> &)>
524     &                         predicate,
525   const VectorTools::NormType norm_type,
526   const double                smallest_abs_coefficient)
527 {
528   Assert(smallest_abs_coefficient >= 0.,
529          ExcMessage("smallest_abs_coefficient should be non-negative."));
530 
531   std::vector<unsigned int> predicate_values;
532   std::vector<double>       norm_values;
533 
534   // first, parse all table elements into a map of predicate values and
535   // coefficients. We could have stored (predicate values ->TableIndicies) map,
536   // but its processing would have been much harder later on.
537   std::map<unsigned int, std::vector<CoefficientType>> pred_to_values;
538   internal::FESeriesImplementation::fill_map(coefficients,
539                                              predicate,
540                                              pred_to_values);
541 
542   // now go through the map and populate the @p norm_values based on @p norm:
543   for (const auto &pred_to_value : pred_to_values)
544     {
545       Vector<CoefficientType> values(pred_to_value.second.cbegin(),
546                                      pred_to_value.second.cend());
547 
548       double norm_value = 0;
549       switch (norm_type)
550         {
551           case VectorTools::L2_norm:
552             {
553               norm_value = values.l2_norm();
554               break;
555             }
556           case VectorTools::L1_norm:
557             {
558               norm_value = values.l1_norm();
559               break;
560             }
561           case VectorTools::Linfty_norm:
562             {
563               norm_value = values.linfty_norm();
564               break;
565             }
566           case VectorTools::mean:
567             {
568               norm_value = internal::FESeriesImplementation::complex_mean_value(
569                 values.mean_value());
570               break;
571             }
572           default:
573             AssertThrow(false, ExcNotImplemented());
574             break;
575         }
576 
577       // will use all non-zero coefficients
578       if (std::abs(norm_value) > smallest_abs_coefficient)
579         {
580           predicate_values.push_back(pred_to_value.first);
581           norm_values.push_back(norm_value);
582         }
583     }
584 
585   return std::make_pair(predicate_values, norm_values);
586 }
587 
588 
589 
590 template <int dim, int spacedim>
591 template <class Archive>
592 inline void
save_transformation_matrices(Archive & ar,const unsigned int)593 FESeries::Fourier<dim, spacedim>::save_transformation_matrices(
594   Archive &ar,
595   const unsigned int /*version*/)
596 {
597   // Store information about those resources which have been used to generate
598   // the transformation matrices.
599   // mode vector
600   ar &n_coefficients_per_direction;
601 
602   // finite element collection
603   unsigned int size = fe_collection->size();
604   ar &         size;
605   for (unsigned int i = 0; i < size; ++i)
606     ar &(*fe_collection)[i].get_name();
607 
608   // quadrature collection
609   size = q_collection.size();
610   ar &size;
611   for (unsigned int i = 0; i < size; ++i)
612     ar &q_collection[i];
613 
614   // Store the actual transform matrices.
615   ar &fourier_transform_matrices;
616 }
617 
618 
619 
620 template <int dim, int spacedim>
621 template <class Archive>
622 inline void
load_transformation_matrices(Archive & ar,const unsigned int)623 FESeries::Fourier<dim, spacedim>::load_transformation_matrices(
624   Archive &ar,
625   const unsigned int /*version*/)
626 {
627   // Check whether the currently registered resources are compatible with
628   // the transformation matrices to load.
629   // mode vector
630   std::vector<unsigned int> compare_coefficients;
631   ar &                      compare_coefficients;
632   Assert(compare_coefficients == n_coefficients_per_direction,
633          ExcMessage("A different number of coefficients vector has been used "
634                     "to generate the transformation matrices you are about "
635                     "to load!"));
636 
637   // finite element collection
638   unsigned int size;
639   ar &         size;
640   AssertDimension(size, fe_collection->size());
641   std::string name;
642   for (unsigned int i = 0; i < size; ++i)
643     {
644       ar &name;
645       Assert(name.compare((*fe_collection)[i].get_name()) == 0,
646              ExcMessage("A different FECollection has been used to generate "
647                         "the transformation matrices you are about to load!"));
648     }
649 
650   // quadrature collection
651   ar &size;
652   AssertDimension(size, q_collection.size());
653   Quadrature<dim> quadrature;
654   for (unsigned int i = 0; i < size; ++i)
655     {
656       ar &quadrature;
657       Assert(quadrature == q_collection[i],
658              ExcMessage("A different QCollection has been used to generate "
659                         "the transformation matrices you are about to load!"));
660     }
661 
662   // Restore the transform matrices since all prerequisites are fulfilled.
663   ar &fourier_transform_matrices;
664 }
665 
666 
667 
668 template <int dim, int spacedim>
669 template <class Archive>
670 inline void
save_transformation_matrices(Archive & ar,const unsigned int)671 FESeries::Legendre<dim, spacedim>::save_transformation_matrices(
672   Archive &ar,
673   const unsigned int /*version*/)
674 {
675   // Store information about those resources which have been used to generate
676   // the transformation matrices.
677   // mode vector
678   ar &n_coefficients_per_direction;
679 
680   // finite element collection
681   unsigned int size = fe_collection->size();
682   ar &         size;
683   for (unsigned int i = 0; i < size; ++i)
684     ar &(*fe_collection)[i].get_name();
685 
686   // quadrature collection
687   size = q_collection.size();
688   ar &size;
689   for (unsigned int i = 0; i < size; ++i)
690     ar &q_collection[i];
691 
692   // Store the actual transform matrices.
693   ar &legendre_transform_matrices;
694 }
695 
696 
697 
698 template <int dim, int spacedim>
699 template <class Archive>
700 inline void
load_transformation_matrices(Archive & ar,const unsigned int)701 FESeries::Legendre<dim, spacedim>::load_transformation_matrices(
702   Archive &ar,
703   const unsigned int /*version*/)
704 {
705   // Check whether the currently registered resources are compatible with
706   // the transformation matrices to load.
707   // mode vector
708   std::vector<unsigned int> compare_coefficients;
709   ar &                      compare_coefficients;
710   Assert(compare_coefficients == n_coefficients_per_direction,
711          ExcMessage("A different number of coefficients vector has been used "
712                     "to generate the transformation matrices you are about "
713                     "to load!"));
714 
715   // finite element collection
716   unsigned int size;
717   ar &         size;
718   AssertDimension(size, fe_collection->size());
719   std::string name;
720   for (unsigned int i = 0; i < size; ++i)
721     {
722       ar &name;
723       Assert(name.compare((*fe_collection)[i].get_name()) == 0,
724              ExcMessage("A different FECollection has been used to generate "
725                         "the transformation matrices you are about to load!"));
726     }
727 
728   // quadrature collection
729   ar &size;
730   AssertDimension(size, q_collection.size());
731   Quadrature<dim> quadrature;
732   for (unsigned int i = 0; i < size; ++i)
733     {
734       ar &quadrature;
735       Assert(quadrature == q_collection[i],
736              ExcMessage("A different QCollection has been used to generate "
737                         "the transformation matrices you are about to load!"));
738     }
739 
740   // Restore the transform matrices since all prerequisites are fulfilled.
741   ar &legendre_transform_matrices;
742 }
743 
744 
745 #endif // DOXYGEN
746 
747 DEAL_II_NAMESPACE_CLOSE
748 
749 #endif // dealii_fe_series_h
750