Reference documentation for deal.II version 9.5.0
\(\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\}}\)
Loading...
Searching...
No Matches
fe_series.h
Go to the documentation of this file.
1// ---------------------------------------------------------------------
2//
3// Copyright (C) 2016 - 2023 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
57namespace FESeries
58{
89 template <int dim, int spacedim = dim>
90 class Fourier : public Subscriptor
91 {
92 public:
93 using CoefficientType = typename std::complex<double>;
94
113 Fourier(const std::vector<unsigned int> & n_coefficients_per_direction,
116 const unsigned int component = numbers::invalid_unsigned_int);
117
123 template <typename Number>
124 void
125 calculate(const ::Vector<Number> &local_dof_values,
126 const unsigned int cell_active_fe_index,
127 Table<dim, CoefficientType> & fourier_coefficients);
128
133 unsigned int
134 get_n_coefficients_per_direction(const unsigned int index) const;
135
146 void
148
158 template <class Archive>
159 void
160 save_transformation_matrices(Archive &ar, const unsigned int version);
161
166 template <class Archive>
167 void
168 load_transformation_matrices(Archive &ar, const unsigned int version);
169
173 bool
174 operator==(const Fourier<dim, spacedim> &fourier) const;
175
176 private:
181 const std::vector<unsigned int> n_coefficients_per_direction;
182
187
192
197
201 std::vector<FullMatrix<CoefficientType>> fourier_transform_matrices;
202
206 std::vector<CoefficientType> unrolled_coefficients;
207
212 const unsigned int component;
213 };
214
215
216
259 template <int dim, int spacedim = dim>
260 class Legendre : public Subscriptor
261 {
262 public:
263 using CoefficientType = double;
264
283 Legendre(const std::vector<unsigned int> &n_coefficients_per_direction,
286 const unsigned int component = numbers::invalid_unsigned_int);
287
293 template <typename Number>
294 void
295 calculate(const ::Vector<Number> &local_dof_values,
296 const unsigned int cell_active_fe_index,
297 Table<dim, CoefficientType> & legendre_coefficients);
298
303 unsigned int
304 get_n_coefficients_per_direction(const unsigned int index) const;
305
316 void
318
328 template <class Archive>
329 void
330 save_transformation_matrices(Archive &ar, const unsigned int version);
331
336 template <class Archive>
337 void
338 load_transformation_matrices(Archive &ar, const unsigned int version);
339
343 bool
344 operator==(const Legendre<dim, spacedim> &legendre) const;
345
346 private:
351 const std::vector<unsigned int> n_coefficients_per_direction;
352
357
362
366 std::vector<FullMatrix<CoefficientType>> legendre_transform_matrices;
367
371 std::vector<CoefficientType> unrolled_coefficients;
372
377 const unsigned int component;
378 };
379
380
381
399 template <int dim, typename CoefficientType>
400 std::pair<std::vector<unsigned int>, std::vector<double>>
402 const std::function<std::pair<bool, unsigned int>(
403 const TableIndices<dim> &)> & predicate,
404 const VectorTools::NormType norm_type,
405 const double smallest_abs_coefficient = 1e-10);
406
412 std::pair<double, double>
413 linear_regression(const std::vector<double> &x, const std::vector<double> &y);
414
415} // namespace FESeries
416
421#ifndef DOXYGEN
422
423// ------------------- inline and template functions ----------------
424
425namespace internal
426{
427 namespace FESeriesImplementation
428 {
429 template <int dim, typename CoefficientType>
430 void
431 fill_map_index(
432 const Table<dim, CoefficientType> &coefficients,
433 const TableIndices<dim> & ind,
434 const std::function<
435 std::pair<bool, unsigned int>(const TableIndices<dim> &)> &predicate,
436 std::map<unsigned int, std::vector<CoefficientType>> &pred_to_values)
437 {
438 const std::pair<bool, unsigned int> pred_pair = predicate(ind);
439 // don't add a value if predicate is false
440 if (pred_pair.first == false)
441 return;
442
443 const unsigned int pred_value = pred_pair.second;
444 const CoefficientType &coeff_value = coefficients(ind);
445 // If pred_value is not in the pred_to_values map, the element will be
446 // created. Otherwise a reference to the existing element is returned.
447 pred_to_values[pred_value].push_back(coeff_value);
448 }
449
450
451
452 template <typename CoefficientType>
453 void
454 fill_map(
455 const Table<1, CoefficientType> &coefficients,
456 const std::function<
457 std::pair<bool, unsigned int>(const TableIndices<1> &)> &predicate,
458 std::map<unsigned int, std::vector<CoefficientType>> & pred_to_values)
459 {
460 for (unsigned int i = 0; i < coefficients.size(0); ++i)
461 {
462 const TableIndices<1> ind(i);
463 fill_map_index(coefficients, ind, predicate, pred_to_values);
464 }
465 }
466
467
468
469 template <typename CoefficientType>
470 void
471 fill_map(
472 const Table<2, CoefficientType> &coefficients,
473 const std::function<
474 std::pair<bool, unsigned int>(const TableIndices<2> &)> &predicate,
475 std::map<unsigned int, std::vector<CoefficientType>> & pred_to_values)
476 {
477 for (unsigned int i = 0; i < coefficients.size(0); ++i)
478 for (unsigned int j = 0; j < coefficients.size(1); ++j)
479 {
480 const TableIndices<2> ind(i, j);
481 fill_map_index(coefficients, ind, predicate, pred_to_values);
482 }
483 }
484
485
486
487 template <typename CoefficientType>
488 void
489 fill_map(
490 const Table<3, CoefficientType> &coefficients,
491 const std::function<
492 std::pair<bool, unsigned int>(const TableIndices<3> &)> &predicate,
493 std::map<unsigned int, std::vector<CoefficientType>> & pred_to_values)
494 {
495 for (unsigned int i = 0; i < coefficients.size(0); ++i)
496 for (unsigned int j = 0; j < coefficients.size(1); ++j)
497 for (unsigned int k = 0; k < coefficients.size(2); ++k)
498 {
499 const TableIndices<3> ind(i, j, k);
500 fill_map_index(coefficients, ind, predicate, pred_to_values);
501 }
502 }
503
504
505
506 template <typename Number>
507 double
508 complex_mean_value(const Number &value)
509 {
510 return value;
511 }
512
513
514
515 template <typename Number>
516 double
517 complex_mean_value(const std::complex<Number> &value)
518 {
519 AssertThrow(false,
521 "FESeries::process_coefficients() can not be used with "
522 "complex-valued coefficients and VectorTools::mean norm."));
523 return std::abs(value);
524 }
525 } // namespace FESeriesImplementation
526} // namespace internal
527
528
529
530template <int dim, typename CoefficientType>
531std::pair<std::vector<unsigned int>, std::vector<double>>
533 const Table<dim, CoefficientType> &coefficients,
534 const std::function<std::pair<bool, unsigned int>(const TableIndices<dim> &)>
535 & predicate,
536 const VectorTools::NormType norm_type,
537 const double smallest_abs_coefficient)
538{
539 Assert(smallest_abs_coefficient >= 0.,
540 ExcMessage("smallest_abs_coefficient should be non-negative."));
541
542 std::vector<unsigned int> predicate_values;
543 std::vector<double> norm_values;
544
545 // first, parse all table elements into a map of predicate values and
546 // coefficients. We could have stored (predicate values ->TableIndices) map,
547 // but its processing would have been much harder later on.
548 std::map<unsigned int, std::vector<CoefficientType>> pred_to_values;
549 internal::FESeriesImplementation::fill_map(coefficients,
550 predicate,
551 pred_to_values);
552
553 // now go through the map and populate the @p norm_values based on @p norm:
554 for (const auto &pred_to_value : pred_to_values)
555 {
556 Vector<CoefficientType> values(pred_to_value.second.cbegin(),
557 pred_to_value.second.cend());
558
559 double norm_value = 0;
560 switch (norm_type)
561 {
563 {
564 norm_value = values.l2_norm();
565 break;
566 }
568 {
569 norm_value = values.l1_norm();
570 break;
571 }
573 {
574 norm_value = values.linfty_norm();
575 break;
576 }
578 {
579 norm_value = internal::FESeriesImplementation::complex_mean_value(
580 values.mean_value());
581 break;
582 }
583 default:
585 break;
586 }
587
588 // will use all non-zero coefficients
589 if (std::abs(norm_value) > smallest_abs_coefficient)
590 {
591 predicate_values.push_back(pred_to_value.first);
592 norm_values.push_back(norm_value);
593 }
594 }
595
596 return std::make_pair(predicate_values, norm_values);
597}
598
599
600
601template <int dim, int spacedim>
602template <class Archive>
603inline void
605 Archive &ar,
606 const unsigned int /*version*/)
607{
608 // Store information about those resources which have been used to generate
609 // the transformation matrices.
610 // mode vector
611 ar &n_coefficients_per_direction;
612
613 // finite element collection
614 unsigned int size = fe_collection->size();
615 ar & size;
616 for (unsigned int i = 0; i < size; ++i)
617 ar &(*fe_collection)[i].get_name();
618
619 // quadrature collection
620 size = q_collection.size();
621 ar &size;
622 for (unsigned int i = 0; i < size; ++i)
623 ar &q_collection[i];
624
625 // Store the actual transform matrices.
626 ar &fourier_transform_matrices;
627}
628
629
630
631template <int dim, int spacedim>
632template <class Archive>
633inline void
635 Archive &ar,
636 const unsigned int /*version*/)
637{
638 // Check whether the currently registered resources are compatible with
639 // the transformation matrices to load.
640 // mode vector
641 std::vector<unsigned int> compare_coefficients;
642 ar & compare_coefficients;
643 Assert(compare_coefficients == n_coefficients_per_direction,
644 ExcMessage("A different number of coefficients vector has been used "
645 "to generate the transformation matrices you are about "
646 "to load!"));
647
648 // finite element collection
649 unsigned int size;
650 ar & size;
651 AssertDimension(size, fe_collection->size());
652 std::string name;
653 for (unsigned int i = 0; i < size; ++i)
654 {
655 ar &name;
656 Assert(name.compare((*fe_collection)[i].get_name()) == 0,
657 ExcMessage("A different FECollection has been used to generate "
658 "the transformation matrices you are about to load!"));
659 }
660
661 // quadrature collection
662 ar &size;
663 AssertDimension(size, q_collection.size());
664 Quadrature<dim> quadrature;
665 for (unsigned int i = 0; i < size; ++i)
666 {
667 ar &quadrature;
668 Assert(quadrature == q_collection[i],
669 ExcMessage("A different QCollection has been used to generate "
670 "the transformation matrices you are about to load!"));
671 }
672
673 // Restore the transform matrices since all prerequisites are fulfilled.
674 ar &fourier_transform_matrices;
675}
676
677
678
679template <int dim, int spacedim>
680template <class Archive>
681inline void
683 Archive &ar,
684 const unsigned int /*version*/)
685{
686 // Store information about those resources which have been used to generate
687 // the transformation matrices.
688 // mode vector
689 ar &n_coefficients_per_direction;
690
691 // finite element collection
692 unsigned int size = fe_collection->size();
693 ar & size;
694 for (unsigned int i = 0; i < size; ++i)
695 ar &(*fe_collection)[i].get_name();
696
697 // quadrature collection
698 size = q_collection.size();
699 ar &size;
700 for (unsigned int i = 0; i < size; ++i)
701 ar &q_collection[i];
702
703 // Store the actual transform matrices.
704 ar &legendre_transform_matrices;
705}
706
707
708
709template <int dim, int spacedim>
710template <class Archive>
711inline void
713 Archive &ar,
714 const unsigned int /*version*/)
715{
716 // Check whether the currently registered resources are compatible with
717 // the transformation matrices to load.
718 // mode vector
719 std::vector<unsigned int> compare_coefficients;
720 ar & compare_coefficients;
721 Assert(compare_coefficients == n_coefficients_per_direction,
722 ExcMessage("A different number of coefficients vector has been used "
723 "to generate the transformation matrices you are about "
724 "to load!"));
725
726 // finite element collection
727 unsigned int size;
728 ar & size;
729 AssertDimension(size, fe_collection->size());
730 std::string name;
731 for (unsigned int i = 0; i < size; ++i)
732 {
733 ar &name;
734 Assert(name.compare((*fe_collection)[i].get_name()) == 0,
735 ExcMessage("A different FECollection has been used to generate "
736 "the transformation matrices you are about to load!"));
737 }
738
739 // quadrature collection
740 ar &size;
741 AssertDimension(size, q_collection.size());
742 Quadrature<dim> quadrature;
743 for (unsigned int i = 0; i < size; ++i)
744 {
745 ar &quadrature;
746 Assert(quadrature == q_collection[i],
747 ExcMessage("A different QCollection has been used to generate "
748 "the transformation matrices you are about to load!"));
749 }
750
751 // Restore the transform matrices since all prerequisites are fulfilled.
752 ar &legendre_transform_matrices;
753}
754
755
756#endif // DOXYGEN
757
759
760#endif // dealii_fe_series_h
SmartPointer< const hp::FECollection< dim, spacedim > > fe_collection
Definition fe_series.h:186
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:196
std::vector< CoefficientType > unrolled_coefficients
Definition fe_series.h:206
const unsigned int component
Definition fe_series.h:212
typename std::complex< double > CoefficientType
Definition fe_series.h:93
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:181
bool operator==(const Fourier< dim, spacedim > &fourier) const
std::vector< FullMatrix< CoefficientType > > fourier_transform_matrices
Definition fe_series.h:201
const hp::QCollection< dim > q_collection
Definition fe_series.h:191
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:351
void save_transformation_matrices(Archive &ar, const unsigned int version)
bool operator==(const Legendre< dim, spacedim > &legendre) const
const unsigned int component
Definition fe_series.h:377
std::vector< FullMatrix< CoefficientType > > legendre_transform_matrices
Definition fe_series.h:366
const hp::QCollection< dim > q_collection
Definition fe_series.h:361
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:371
SmartPointer< const hp::FECollection< dim, spacedim > > fe_collection
Definition fe_series.h:356
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)
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:472
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:473
static ::ExceptionBase & ExcNotImplemented()
#define Assert(cond, exc)
#define AssertDimension(dim1, dim2)
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
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:213
::VectorizedArray< Number, width > abs(const ::VectorizedArray< Number, width > &)