deal.II version GIT relicensing-1982-gbb401394a0 2024-10-13 08:40:01+00:00
\(\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// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 2016 - 2023 by the deal.II authors
5//
6// This file is part of the deal.II library.
7//
8// Part of the source code is dual licensed under Apache-2.0 WITH
9// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
10// governing the source code and code contributions can be found in
11// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
12//
13// ------------------------------------------------------------------------
14
15#ifndef dealii_fe_series_h
16#define dealii_fe_series_h
17
18
19
20#include <deal.II/base/config.h>
21
24#include <deal.II/base/table.h>
26#include <deal.II/base/tensor.h>
27
30
32#include <deal.II/lac/vector.h>
33
35
36#include <memory>
37#include <string>
38#include <vector>
39
40
42
43
56namespace FESeries
57{
88 template <int dim, int spacedim = dim>
90 {
91 public:
92 using CoefficientType = typename std::complex<double>;
93
112 Fourier(const std::vector<unsigned int> &n_coefficients_per_direction,
115 const unsigned int component = numbers::invalid_unsigned_int);
116
122 template <typename Number>
123 void
124 calculate(const ::Vector<Number> &local_dof_values,
125 const unsigned int cell_active_fe_index,
126 Table<dim, CoefficientType> &fourier_coefficients);
127
132 unsigned int
133 get_n_coefficients_per_direction(const unsigned int index) const;
134
145 void
147
157 template <class Archive>
158 void
159 save_transformation_matrices(Archive &ar, const unsigned int version);
160
165 template <class Archive>
166 void
167 load_transformation_matrices(Archive &ar, const unsigned int version);
168
172 bool
173 operator==(const Fourier<dim, spacedim> &fourier) const;
174
175 private:
180 const std::vector<unsigned int> n_coefficients_per_direction;
181
186
191
196
200 std::vector<FullMatrix<CoefficientType>> fourier_transform_matrices;
201
205 std::vector<CoefficientType> unrolled_coefficients;
206
211 const unsigned int component;
212 };
213
214
215
258 template <int dim, int spacedim = dim>
260 {
261 public:
262 using CoefficientType = double;
263
282 Legendre(const std::vector<unsigned int> &n_coefficients_per_direction,
285 const unsigned int component = numbers::invalid_unsigned_int);
286
292 template <typename Number>
293 void
294 calculate(const ::Vector<Number> &local_dof_values,
295 const unsigned int cell_active_fe_index,
296 Table<dim, CoefficientType> &legendre_coefficients);
297
302 unsigned int
303 get_n_coefficients_per_direction(const unsigned int index) const;
304
315 void
317
327 template <class Archive>
328 void
329 save_transformation_matrices(Archive &ar, const unsigned int version);
330
335 template <class Archive>
336 void
337 load_transformation_matrices(Archive &ar, const unsigned int version);
338
342 bool
343 operator==(const Legendre<dim, spacedim> &legendre) const;
344
345 private:
350 const std::vector<unsigned int> n_coefficients_per_direction;
351
356
361
365 std::vector<FullMatrix<CoefficientType>> legendre_transform_matrices;
366
370 std::vector<CoefficientType> unrolled_coefficients;
371
376 const unsigned int component;
377 };
378
379
380
398 template <int dim, typename CoefficientType>
399 std::pair<std::vector<unsigned int>, std::vector<double>>
401 const std::function<std::pair<bool, unsigned int>(
402 const TableIndices<dim> &)> &predicate,
403 const VectorTools::NormType norm_type,
404 const double smallest_abs_coefficient = 1e-10);
405
411 std::pair<double, double>
412 linear_regression(const std::vector<double> &x, const std::vector<double> &y);
413
414} // namespace FESeries
415
420#ifndef DOXYGEN
421
422// ------------------- inline and template functions ----------------
423
424namespace internal
425{
426 namespace FESeriesImplementation
427 {
428 template <int dim, typename CoefficientType>
429 void
430 fill_map_index(
431 const Table<dim, CoefficientType> &coefficients,
432 const TableIndices<dim> &ind,
433 const std::function<
434 std::pair<bool, unsigned int>(const TableIndices<dim> &)> &predicate,
435 std::map<unsigned int, std::vector<CoefficientType>> &pred_to_values)
436 {
437 const std::pair<bool, unsigned int> pred_pair = predicate(ind);
438 // don't add a value if predicate is false
439 if (pred_pair.first == false)
440 return;
441
442 const unsigned int pred_value = pred_pair.second;
443 const CoefficientType &coeff_value = coefficients(ind);
444 // If pred_value is not in the pred_to_values map, the element will be
445 // created. Otherwise a reference to the existing element is returned.
446 pred_to_values[pred_value].push_back(coeff_value);
447 }
448
449
450
451 template <typename CoefficientType>
452 void
453 fill_map(
454 const Table<1, CoefficientType> &coefficients,
455 const std::function<
456 std::pair<bool, unsigned int>(const TableIndices<1> &)> &predicate,
457 std::map<unsigned int, std::vector<CoefficientType>> &pred_to_values)
458 {
459 for (unsigned int i = 0; i < coefficients.size(0); ++i)
460 {
461 const TableIndices<1> ind(i);
462 fill_map_index(coefficients, ind, predicate, pred_to_values);
463 }
464 }
465
466
467
468 template <typename CoefficientType>
469 void
470 fill_map(
471 const Table<2, CoefficientType> &coefficients,
472 const std::function<
473 std::pair<bool, unsigned int>(const TableIndices<2> &)> &predicate,
474 std::map<unsigned int, std::vector<CoefficientType>> &pred_to_values)
475 {
476 for (unsigned int i = 0; i < coefficients.size(0); ++i)
477 for (unsigned int j = 0; j < coefficients.size(1); ++j)
478 {
479 const TableIndices<2> ind(i, j);
480 fill_map_index(coefficients, ind, predicate, pred_to_values);
481 }
482 }
483
484
485
486 template <typename CoefficientType>
487 void
488 fill_map(
489 const Table<3, CoefficientType> &coefficients,
490 const std::function<
491 std::pair<bool, unsigned int>(const TableIndices<3> &)> &predicate,
492 std::map<unsigned int, std::vector<CoefficientType>> &pred_to_values)
493 {
494 for (unsigned int i = 0; i < coefficients.size(0); ++i)
495 for (unsigned int j = 0; j < coefficients.size(1); ++j)
496 for (unsigned int k = 0; k < coefficients.size(2); ++k)
497 {
498 const TableIndices<3> ind(i, j, k);
499 fill_map_index(coefficients, ind, predicate, pred_to_values);
500 }
501 }
502
503
504
505 template <typename Number>
506 double
507 complex_mean_value(const Number &value)
508 {
509 return value;
510 }
511
512
513
514 template <typename Number>
515 double
516 complex_mean_value(const std::complex<Number> &value)
517 {
518 AssertThrow(false,
520 "FESeries::process_coefficients() can not be used with "
521 "complex-valued coefficients and VectorTools::mean norm."));
522 return std::abs(value);
523 }
524 } // namespace FESeriesImplementation
525} // namespace internal
526
527
528
529template <int dim, typename CoefficientType>
530std::pair<std::vector<unsigned int>, std::vector<double>>
532 const Table<dim, CoefficientType> &coefficients,
533 const std::function<std::pair<bool, unsigned int>(const TableIndices<dim> &)>
534 &predicate,
535 const VectorTools::NormType norm_type,
536 const double smallest_abs_coefficient)
537{
538 Assert(smallest_abs_coefficient >= 0.,
539 ExcMessage("smallest_abs_coefficient should be non-negative."));
540
541 std::vector<unsigned int> predicate_values;
542 std::vector<double> norm_values;
543
544 // first, parse all table elements into a map of predicate values and
545 // coefficients. We could have stored (predicate values ->TableIndices) map,
546 // but its processing would have been much harder later on.
547 std::map<unsigned int, std::vector<CoefficientType>> pred_to_values;
548 internal::FESeriesImplementation::fill_map(coefficients,
549 predicate,
550 pred_to_values);
551
552 // now go through the map and populate the @p norm_values based on @p norm:
553 for (const auto &pred_to_value : pred_to_values)
554 {
555 Vector<CoefficientType> values(pred_to_value.second.cbegin(),
556 pred_to_value.second.cend());
557
558 double norm_value = 0;
559 switch (norm_type)
560 {
562 {
563 norm_value = values.l2_norm();
564 break;
565 }
567 {
568 norm_value = values.l1_norm();
569 break;
570 }
572 {
573 norm_value = values.linfty_norm();
574 break;
575 }
577 {
578 norm_value = internal::FESeriesImplementation::complex_mean_value(
579 values.mean_value());
580 break;
581 }
582 default:
584 break;
585 }
586
587 // will use all non-zero coefficients
588 if (std::abs(norm_value) > smallest_abs_coefficient)
589 {
590 predicate_values.push_back(pred_to_value.first);
591 norm_values.push_back(norm_value);
592 }
593 }
594
595 return std::make_pair(predicate_values, norm_values);
596}
597
598
599
600template <int dim, int spacedim>
601template <class Archive>
602inline void
604 Archive &ar,
605 const unsigned int /*version*/)
606{
607 // Store information about those resources which have been used to generate
608 // the transformation matrices.
609 // mode vector
610 ar &n_coefficients_per_direction;
611
612 // finite element collection
613 unsigned int size = fe_collection->size();
614 ar &size;
615 for (unsigned int i = 0; i < size; ++i)
616 ar &(*fe_collection)[i].get_name();
617
618 // quadrature collection
619 size = q_collection.size();
620 ar &size;
621 for (unsigned int i = 0; i < size; ++i)
622 ar &q_collection[i];
623
624 // Store the actual transform matrices.
625 ar &fourier_transform_matrices;
626}
627
628
629
630template <int dim, int spacedim>
631template <class Archive>
632inline void
634 Archive &ar,
635 const unsigned int /*version*/)
636{
637 // Check whether the currently registered resources are compatible with
638 // the transformation matrices to load.
639 // mode vector
640 std::vector<unsigned int> compare_coefficients;
641 ar &compare_coefficients;
642 Assert(compare_coefficients == n_coefficients_per_direction,
643 ExcMessage("A different number of coefficients vector has been used "
644 "to generate the transformation matrices you are about "
645 "to load!"));
646
647 // finite element collection
648 unsigned int size;
649 ar &size;
650 AssertDimension(size, fe_collection->size());
651 std::string name;
652 for (unsigned int i = 0; i < size; ++i)
653 {
654 ar &name;
655 Assert(name.compare((*fe_collection)[i].get_name()) == 0,
656 ExcMessage("A different FECollection has been used to generate "
657 "the transformation matrices you are about to load!"));
658 }
659
660 // quadrature collection
661 ar &size;
662 AssertDimension(size, q_collection.size());
663 Quadrature<dim> quadrature;
664 for (unsigned int i = 0; i < size; ++i)
665 {
666 ar &quadrature;
667 Assert(quadrature == q_collection[i],
668 ExcMessage("A different QCollection has been used to generate "
669 "the transformation matrices you are about to load!"));
670 }
671
672 // Restore the transform matrices since all prerequisites are fulfilled.
673 ar &fourier_transform_matrices;
674}
675
676
677
678template <int dim, int spacedim>
679template <class Archive>
680inline void
682 Archive &ar,
683 const unsigned int /*version*/)
684{
685 // Store information about those resources which have been used to generate
686 // the transformation matrices.
687 // mode vector
688 ar &n_coefficients_per_direction;
689
690 // finite element collection
691 unsigned int size = fe_collection->size();
692 ar &size;
693 for (unsigned int i = 0; i < size; ++i)
694 ar &(*fe_collection)[i].get_name();
695
696 // quadrature collection
697 size = q_collection.size();
698 ar &size;
699 for (unsigned int i = 0; i < size; ++i)
700 ar &q_collection[i];
701
702 // Store the actual transform matrices.
703 ar &legendre_transform_matrices;
704}
705
706
707
708template <int dim, int spacedim>
709template <class Archive>
710inline void
712 Archive &ar,
713 const unsigned int /*version*/)
714{
715 // Check whether the currently registered resources are compatible with
716 // the transformation matrices to load.
717 // mode vector
718 std::vector<unsigned int> compare_coefficients;
719 ar &compare_coefficients;
720 Assert(compare_coefficients == n_coefficients_per_direction,
721 ExcMessage("A different number of coefficients vector has been used "
722 "to generate the transformation matrices you are about "
723 "to load!"));
724
725 // finite element collection
726 unsigned int size;
727 ar &size;
728 AssertDimension(size, fe_collection->size());
729 std::string name;
730 for (unsigned int i = 0; i < size; ++i)
731 {
732 ar &name;
733 Assert(name.compare((*fe_collection)[i].get_name()) == 0,
734 ExcMessage("A different FECollection has been used to generate "
735 "the transformation matrices you are about to load!"));
736 }
737
738 // quadrature collection
739 ar &size;
740 AssertDimension(size, q_collection.size());
741 Quadrature<dim> quadrature;
742 for (unsigned int i = 0; i < size; ++i)
743 {
744 ar &quadrature;
745 Assert(quadrature == q_collection[i],
746 ExcMessage("A different QCollection has been used to generate "
747 "the transformation matrices you are about to load!"));
748 }
749
750 // Restore the transform matrices since all prerequisites are fulfilled.
751 ar &legendre_transform_matrices;
752}
753
754
755#endif // DOXYGEN
756
758
759#endif // dealii_fe_series_h
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:195
std::vector< CoefficientType > unrolled_coefficients
Definition fe_series.h:205
const unsigned int component
Definition fe_series.h:211
typename std::complex< double > CoefficientType
Definition fe_series.h:92
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:180
ObserverPointer< const hp::FECollection< dim, spacedim > > fe_collection
Definition fe_series.h:185
bool operator==(const Fourier< dim, spacedim > &fourier) const
std::vector< FullMatrix< CoefficientType > > fourier_transform_matrices
Definition fe_series.h:200
const hp::QCollection< dim > q_collection
Definition fe_series.h:190
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:350
ObserverPointer< const hp::FECollection< dim, spacedim > > fe_collection
Definition fe_series.h:355
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:376
std::vector< FullMatrix< CoefficientType > > legendre_transform_matrices
Definition fe_series.h:365
const hp::QCollection< dim > q_collection
Definition fe_series.h:360
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:370
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:498
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:499
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:29
static const unsigned int invalid_unsigned_int
Definition types.h:220
::VectorizedArray< Number, width > abs(const ::VectorizedArray< Number, width > &)