Loading [MathJax]/extensions/TeX/newcommand.js
 Reference documentation for deal.II version 9.4.1
\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}} \newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=} \newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]} \newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Modules Pages
fe_series.h
Go to the documentation of this file.
1// ---------------------------------------------------------------------
2//
3// Copyright (C) 2016 - 2021 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
25#include <deal.II/base/table.h>
27#include <deal.II/base/tensor.h>
28
31
33#include <deal.II/lac/vector.h>
34
36
37#include <memory>
38#include <string>
39#include <vector>
40
41
43
44
47
48
55namespace FESeries
56{
87 template <int dim, int spacedim = dim>
88 class Fourier : public Subscriptor
89 {
90 public:
91 using CoefficientType = typename std::complex<double>;
92
111 Fourier(const std::vector<unsigned int> & n_coefficients_per_direction,
114 const unsigned int component = numbers::invalid_unsigned_int);
115
121 template <typename Number>
122 void
123 calculate(const ::Vector<Number> &local_dof_values,
124 const unsigned int cell_active_fe_index,
125 Table<dim, CoefficientType> & fourier_coefficients);
126
131 unsigned int
132 get_n_coefficients_per_direction(const unsigned int index) const;
133
144 void
146
156 template <class Archive>
157 void
158 save_transformation_matrices(Archive &ar, const unsigned int version);
159
164 template <class Archive>
165 void
166 load_transformation_matrices(Archive &ar, const unsigned int version);
167
171 bool
172 operator==(const Fourier<dim, spacedim> &fourier) const;
173
174 private:
179 const std::vector<unsigned int> n_coefficients_per_direction;
180
185
190
195
199 std::vector<FullMatrix<CoefficientType>> fourier_transform_matrices;
200
204 std::vector<CoefficientType> unrolled_coefficients;
205
210 const unsigned int component;
211 };
212
213
214
257 template <int dim, int spacedim = dim>
258 class Legendre : public Subscriptor
259 {
260 public:
261 using CoefficientType = double;
262
281 Legendre(const std::vector<unsigned int> &n_coefficients_per_direction,
284 const unsigned int component = numbers::invalid_unsigned_int);
285
291 template <typename Number>
292 void
293 calculate(const ::Vector<Number> &local_dof_values,
294 const unsigned int cell_active_fe_index,
295 Table<dim, CoefficientType> & legendre_coefficients);
296
301 unsigned int
302 get_n_coefficients_per_direction(const unsigned int index) const;
303
314 void
316
326 template <class Archive>
327 void
328 save_transformation_matrices(Archive &ar, const unsigned int version);
329
334 template <class Archive>
335 void
336 load_transformation_matrices(Archive &ar, const unsigned int version);
337
341 bool
342 operator==(const Legendre<dim, spacedim> &legendre) const;
343
344 private:
349 const std::vector<unsigned int> n_coefficients_per_direction;
350
355
360
364 std::vector<FullMatrix<CoefficientType>> legendre_transform_matrices;
365
369 std::vector<CoefficientType> unrolled_coefficients;
370
375 const unsigned int component;
376 };
377
378
379
397 template <int dim, typename CoefficientType>
398 std::pair<std::vector<unsigned int>, std::vector<double>>
400 const std::function<std::pair<bool, unsigned int>(
401 const TableIndices<dim> &)> & predicate,
402 const VectorTools::NormType norm_type,
403 const double smallest_abs_coefficient = 1e-10);
404
410 std::pair<double, double>
411 linear_regression(const std::vector<double> &x, const std::vector<double> &y);
412
413} // namespace FESeries
414
419#ifndef DOXYGEN
420
421// ------------------- inline and template functions ----------------
422
423namespace internal
424{
425 namespace FESeriesImplementation
426 {
427 template <int dim, typename CoefficientType>
428 void
429 fill_map_index(
430 const Table<dim, CoefficientType> &coefficients,
431 const TableIndices<dim> & ind,
432 const std::function<
433 std::pair<bool, unsigned int>(const TableIndices<dim> &)> &predicate,
434 std::map<unsigned int, std::vector<CoefficientType>> &pred_to_values)
435 {
436 const std::pair<bool, unsigned int> pred_pair = predicate(ind);
437 // don't add a value if predicate is false
438 if (pred_pair.first == false)
439 return;
440
441 const unsigned int pred_value = pred_pair.second;
442 const CoefficientType &coeff_value = coefficients(ind);
443 // If pred_value is not in the pred_to_values map, the element will be
444 // created. Otherwise a reference to the existing element is returned.
445 pred_to_values[pred_value].push_back(coeff_value);
446 }
447
448
449
450 template <typename CoefficientType>
451 void
452 fill_map(
453 const Table<1, CoefficientType> &coefficients,
454 const std::function<
455 std::pair<bool, unsigned int>(const TableIndices<1> &)> &predicate,
456 std::map<unsigned int, std::vector<CoefficientType>> & pred_to_values)
457 {
458 for (unsigned int i = 0; i < coefficients.size(0); ++i)
459 {
460 const TableIndices<1> ind(i);
461 fill_map_index(coefficients, ind, predicate, pred_to_values);
462 }
463 }
464
465
466
467 template <typename CoefficientType>
468 void
469 fill_map(
470 const Table<2, CoefficientType> &coefficients,
471 const std::function<
472 std::pair<bool, unsigned int>(const TableIndices<2> &)> &predicate,
473 std::map<unsigned int, std::vector<CoefficientType>> & pred_to_values)
474 {
475 for (unsigned int i = 0; i < coefficients.size(0); ++i)
476 for (unsigned int j = 0; j < coefficients.size(1); ++j)
477 {
478 const TableIndices<2> ind(i, j);
479 fill_map_index(coefficients, ind, predicate, pred_to_values);
480 }
481 }
482
483
484
485 template <typename CoefficientType>
486 void
487 fill_map(
488 const Table<3, CoefficientType> &coefficients,
489 const std::function<
490 std::pair<bool, unsigned int>(const TableIndices<3> &)> &predicate,
491 std::map<unsigned int, std::vector<CoefficientType>> & pred_to_values)
492 {
493 for (unsigned int i = 0; i < coefficients.size(0); ++i)
494 for (unsigned int j = 0; j < coefficients.size(1); ++j)
495 for (unsigned int k = 0; k < coefficients.size(2); ++k)
496 {
497 const TableIndices<3> ind(i, j, k);
498 fill_map_index(coefficients, ind, predicate, pred_to_values);
499 }
500 }
501
502
503
504 template <typename Number>
505 double
506 complex_mean_value(const Number &value)
507 {
508 return value;
509 }
510
511
512
513 template <typename Number>
514 double
515 complex_mean_value(const std::complex<Number> &value)
516 {
517 AssertThrow(false,
519 "FESeries::process_coefficients() can not be used with "
520 "complex-valued coefficients and VectorTools::mean norm."));
521 return std::abs(value);
522 }
523 } // namespace FESeriesImplementation
524} // namespace internal
525
526
527
528template <int dim, typename CoefficientType>
529std::pair<std::vector<unsigned int>, std::vector<double>>
531 const Table<dim, CoefficientType> &coefficients,
532 const std::function<std::pair<bool, unsigned int>(const TableIndices<dim> &)>
533 & predicate,
534 const VectorTools::NormType norm_type,
535 const double smallest_abs_coefficient)
536{
537 Assert(smallest_abs_coefficient >= 0.,
538 ExcMessage("smallest_abs_coefficient should be non-negative."));
539
540 std::vector<unsigned int> predicate_values;
541 std::vector<double> norm_values;
542
543 // first, parse all table elements into a map of predicate values and
544 // coefficients. We could have stored (predicate values ->TableIndicies) map,
545 // but its processing would have been much harder later on.
546 std::map<unsigned int, std::vector<CoefficientType>> pred_to_values;
547 internal::FESeriesImplementation::fill_map(coefficients,
548 predicate,
549 pred_to_values);
550
551 // now go through the map and populate the @p norm_values based on @p norm:
552 for (const auto &pred_to_value : pred_to_values)
553 {
554 Vector<CoefficientType> values(pred_to_value.second.cbegin(),
555 pred_to_value.second.cend());
556
557 double norm_value = 0;
558 switch (norm_type)
559 {
561 {
562 norm_value = values.l2_norm();
563 break;
564 }
566 {
567 norm_value = values.l1_norm();
568 break;
569 }
571 {
572 norm_value = values.linfty_norm();
573 break;
574 }
576 {
577 norm_value = internal::FESeriesImplementation::complex_mean_value(
578 values.mean_value());
579 break;
580 }
581 default:
583 break;
584 }
585
586 // will use all non-zero coefficients
587 if (std::abs(norm_value) > smallest_abs_coefficient)
588 {
589 predicate_values.push_back(pred_to_value.first);
590 norm_values.push_back(norm_value);
591 }
592 }
593
594 return std::make_pair(predicate_values, norm_values);
595}
596
597
598
599template <int dim, int spacedim>
600template <class Archive>
601inline void
603 Archive &ar,
604 const unsigned int /*version*/)
605{
606 // Store information about those resources which have been used to generate
607 // the transformation matrices.
608 // mode vector
609 ar &n_coefficients_per_direction;
610
611 // finite element collection
612 unsigned int size = fe_collection->size();
613 ar & size;
614 for (unsigned int i = 0; i < size; ++i)
615 ar &(*fe_collection)[i].get_name();
616
617 // quadrature collection
618 size = q_collection.size();
619 ar &size;
620 for (unsigned int i = 0; i < size; ++i)
621 ar &q_collection[i];
622
623 // Store the actual transform matrices.
624 ar &fourier_transform_matrices;
625}
626
627
628
629template <int dim, int spacedim>
630template <class Archive>
631inline void
633 Archive &ar,
634 const unsigned int /*version*/)
635{
636 // Check whether the currently registered resources are compatible with
637 // the transformation matrices to load.
638 // mode vector
639 std::vector<unsigned int> compare_coefficients;
640 ar & compare_coefficients;
641 Assert(compare_coefficients == n_coefficients_per_direction,
642 ExcMessage("A different number of coefficients vector has been used "
643 "to generate the transformation matrices you are about "
644 "to load!"));
645
646 // finite element collection
647 unsigned int size;
648 ar & size;
649 AssertDimension(size, fe_collection->size());
650 std::string name;
651 for (unsigned int i = 0; i < size; ++i)
652 {
653 ar &name;
654 Assert(name.compare((*fe_collection)[i].get_name()) == 0,
655 ExcMessage("A different FECollection has been used to generate "
656 "the transformation matrices you are about to load!"));
657 }
658
659 // quadrature collection
660 ar &size;
661 AssertDimension(size, q_collection.size());
662 Quadrature<dim> quadrature;
663 for (unsigned int i = 0; i < size; ++i)
664 {
665 ar &quadrature;
666 Assert(quadrature == q_collection[i],
667 ExcMessage("A different QCollection has been used to generate "
668 "the transformation matrices you are about to load!"));
669 }
670
671 // Restore the transform matrices since all prerequisites are fulfilled.
672 ar &fourier_transform_matrices;
673}
674
675
676
677template <int dim, int spacedim>
678template <class Archive>
679inline void
681 Archive &ar,
682 const unsigned int /*version*/)
683{
684 // Store information about those resources which have been used to generate
685 // the transformation matrices.
686 // mode vector
687 ar &n_coefficients_per_direction;
688
689 // finite element collection
690 unsigned int size = fe_collection->size();
691 ar & size;
692 for (unsigned int i = 0; i < size; ++i)
693 ar &(*fe_collection)[i].get_name();
694
695 // quadrature collection
696 size = q_collection.size();
697 ar &size;
698 for (unsigned int i = 0; i < size; ++i)
699 ar &q_collection[i];
700
701 // Store the actual transform matrices.
702 ar &legendre_transform_matrices;
703}
704
705
706
707template <int dim, int spacedim>
708template <class Archive>
709inline void
711 Archive &ar,
712 const unsigned int /*version*/)
713{
714 // Check whether the currently registered resources are compatible with
715 // the transformation matrices to load.
716 // mode vector
717 std::vector<unsigned int> compare_coefficients;
718 ar & compare_coefficients;
719 Assert(compare_coefficients == n_coefficients_per_direction,
720 ExcMessage("A different number of coefficients vector has been used "
721 "to generate the transformation matrices you are about "
722 "to load!"));
723
724 // finite element collection
725 unsigned int size;
726 ar & size;
727 AssertDimension(size, fe_collection->size());
728 std::string name;
729 for (unsigned int i = 0; i < size; ++i)
730 {
731 ar &name;
732 Assert(name.compare((*fe_collection)[i].get_name()) == 0,
733 ExcMessage("A different FECollection has been used to generate "
734 "the transformation matrices you are about to load!"));
735 }
736
737 // quadrature collection
738 ar &size;
739 AssertDimension(size, q_collection.size());
740 Quadrature<dim> quadrature;
741 for (unsigned int i = 0; i < size; ++i)
742 {
743 ar &quadrature;
744 Assert(quadrature == q_collection[i],
745 ExcMessage("A different QCollection has been used to generate "
746 "the transformation matrices you are about to load!"));
747 }
748
749 // Restore the transform matrices since all prerequisites are fulfilled.
750 ar &legendre_transform_matrices;
751}
752
753
754#endif // DOXYGEN
755
757
758#endif // dealii_fe_series_h
SmartPointer< const hp::FECollection< dim, spacedim > > fe_collection
Definition: fe_series.h:184
void calculate(const ::Vector< Number > &local_dof_values, const unsigned int cell_active_fe_index, Table< dim, CoefficientType > &fourier_coefficients)
Table< dim, Tensor< 1, dim > > k_vectors
Definition: fe_series.h:194
std::vector< CoefficientType > unrolled_coefficients
Definition: fe_series.h:204
const unsigned int component
Definition: fe_series.h:210
typename std::complex< double > CoefficientType
Definition: fe_series.h:91
unsigned int get_n_coefficients_per_direction(const unsigned int index) const
const std::vector< unsigned int > n_coefficients_per_direction
Definition: fe_series.h:179
bool operator==(const Fourier< dim, spacedim > &fourier) const
std::vector< FullMatrix< CoefficientType > > fourier_transform_matrices
Definition: fe_series.h:199
const hp::QCollection< dim > q_collection
Definition: fe_series.h:189
void precalculate_all_transformation_matrices()
void load_transformation_matrices(Archive &ar, const unsigned int version)
void save_transformation_matrices(Archive &ar, const unsigned int version)
const std::vector< unsigned int > n_coefficients_per_direction
Definition: fe_series.h:349
void save_transformation_matrices(Archive &ar, const unsigned int version)
double CoefficientType
Definition: fe_series.h:261
bool operator==(const Legendre< dim, spacedim > &legendre) const
const unsigned int component
Definition: fe_series.h:375
std::vector< FullMatrix< CoefficientType > > legendre_transform_matrices
Definition: fe_series.h:364
const hp::QCollection< dim > q_collection
Definition: fe_series.h:359
void precalculate_all_transformation_matrices()
unsigned int get_n_coefficients_per_direction(const unsigned int index) const
std::vector< CoefficientType > unrolled_coefficients
Definition: fe_series.h:369
SmartPointer< const hp::FECollection< dim, spacedim > > fe_collection
Definition: fe_series.h:354
void load_transformation_matrices(Archive &ar, const unsigned int version)
void calculate(const ::Vector< Number > &local_dof_values, const unsigned int cell_active_fe_index, Table< dim, CoefficientType > &legendre_coefficients)
Definition: vector.h:109
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:442
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:443
static ::ExceptionBase & ExcNotImplemented()
#define Assert(cond, exc)
Definition: exceptions.h:1473
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1667
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
Definition: exceptions.h:1583
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=1e-10)
std::pair< double, double > linear_regression(const std::vector< double > &x, const std::vector< double > &y)
Definition: fe_series.cc:30
static const unsigned int invalid_unsigned_int
Definition: types.h:201
::VectorizedArray< Number, width > abs(const ::VectorizedArray< Number, width > &)