Reference documentation for deal.II version 9.1.1
\(\newcommand{\dealcoloneq}{\mathrel{\vcenter{:}}=}\)
fe_series_fourier.cc
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2016 - 2019 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 
17 
18 #include <deal.II/base/config.h>
19 
20 #include <deal.II/base/numbers.h>
21 
22 #include <deal.II/fe/fe_series.h>
23 
24 #include <iostream>
25 
26 
27 DEAL_II_NAMESPACE_OPEN
28 
29 namespace
30 {
31  void set_k_vectors(Table<1, Tensor<1, 1>> &k_vectors, const unsigned int N)
32  {
33  k_vectors.reinit(TableIndices<1>(N));
34  for (unsigned int i = 0; i < N; ++i)
35  k_vectors(i)[0] = 2. * numbers::PI * i;
36  }
37 
38  void set_k_vectors(Table<2, Tensor<1, 2>> &k_vectors, const unsigned int N)
39  {
40  k_vectors.reinit(TableIndices<2>(N, N));
41  for (unsigned int i = 0; i < N; ++i)
42  for (unsigned int j = 0; j < N; ++j)
43  {
44  k_vectors(i, j)[0] = 2. * numbers::PI * i;
45  k_vectors(i, j)[1] = 2. * numbers::PI * j;
46  }
47  }
48 
49  void set_k_vectors(Table<3, Tensor<1, 3>> &k_vectors, const unsigned int N)
50  {
51  k_vectors.reinit(TableIndices<3>(N, N, N));
52  for (unsigned int i = 0; i < N; ++i)
53  for (unsigned int j = 0; j < N; ++j)
54  for (unsigned int k = 0; k < N; ++k)
55  {
56  k_vectors(i, j, k)[0] = 2. * numbers::PI * i;
57  k_vectors(i, j, k)[1] = 2. * numbers::PI * j;
58  k_vectors(i, j, k)[2] = 2. * numbers::PI * k;
59  }
60  }
61 
62 
63 
64  template <int dim, int spacedim>
65  std::complex<double>
66  integrate(const FiniteElement<dim, spacedim> &fe,
67  const Quadrature<dim> & quadrature,
68  const Tensor<1, dim> & k_vector,
69  const unsigned int j)
70  {
71  std::complex<double> sum = 0;
72  for (unsigned int q = 0; q < quadrature.size(); ++q)
73  {
74  const Point<dim> &x_q = quadrature.point(q);
75  sum += std::exp(std::complex<double>(0, 1) * (k_vector * x_q)) *
76  fe.shape_value(j, x_q) * quadrature.weight(q);
77  }
78  return sum;
79  }
80 
81 
82 
83  /*
84  * Ensure that the transformation matrix for FiniteElement index
85  * @p fe_index is calculated. If not, calculate it.
86  */
87  template <int spacedim>
88  void
89  ensure_existence(
90  const hp::FECollection<1, spacedim> & fe_collection,
91  const hp::QCollection<1> & q_collection,
92  const Table<1, Tensor<1, 1>> & k_vectors,
93  const unsigned int fe,
94  std::vector<FullMatrix<std::complex<double>>> &fourier_transform_matrices)
95  {
96  AssertIndexRange(fe, fe_collection.size());
97 
98  if (fourier_transform_matrices[fe].m() == 0)
99  {
100  fourier_transform_matrices[fe].reinit(k_vectors.n_elements(),
101  fe_collection[fe].dofs_per_cell);
102 
103  for (unsigned int k = 0; k < k_vectors.size(0); ++k)
104  for (unsigned int j = 0; j < fe_collection[fe].dofs_per_cell; ++j)
105  fourier_transform_matrices[fe](k, j) =
106  integrate(fe_collection[fe], q_collection[fe], k_vectors(k), j);
107  }
108  }
109 
110  template <int spacedim>
111  void
112  ensure_existence(
113  const hp::FECollection<2, spacedim> & fe_collection,
114  const hp::QCollection<2> & q_collection,
115  const Table<2, Tensor<1, 2>> & k_vectors,
116  const unsigned int fe,
117  std::vector<FullMatrix<std::complex<double>>> &fourier_transform_matrices)
118  {
119  AssertIndexRange(fe, fe_collection.size());
120 
121  if (fourier_transform_matrices[fe].m() == 0)
122  {
123  fourier_transform_matrices[fe].reinit(k_vectors.n_elements(),
124  fe_collection[fe].dofs_per_cell);
125 
126  unsigned int k = 0;
127  for (unsigned int k1 = 0; k1 < k_vectors.size(0); ++k1)
128  for (unsigned int k2 = 0; k2 < k_vectors.size(1); ++k2, k++)
129  for (unsigned int j = 0; j < fe_collection[fe].dofs_per_cell; ++j)
130  fourier_transform_matrices[fe](k, j) = integrate(
131  fe_collection[fe], q_collection[fe], k_vectors(k1, k2), j);
132  }
133  }
134 
135  template <int spacedim>
136  void
137  ensure_existence(
138  const hp::FECollection<3, spacedim> & fe_collection,
139  const hp::QCollection<3> & q_collection,
140  const Table<3, Tensor<1, 3>> & k_vectors,
141  const unsigned int fe,
142  std::vector<FullMatrix<std::complex<double>>> &fourier_transform_matrices)
143  {
144  AssertIndexRange(fe, fe_collection.size());
145 
146  if (fourier_transform_matrices[fe].m() == 0)
147  {
148  fourier_transform_matrices[fe].reinit(k_vectors.n_elements(),
149  fe_collection[fe].dofs_per_cell);
150 
151  unsigned int k = 0;
152  for (unsigned int k1 = 0; k1 < k_vectors.size(0); ++k1)
153  for (unsigned int k2 = 0; k2 < k_vectors.size(1); ++k2)
154  for (unsigned int k3 = 0; k3 < k_vectors.size(2); ++k3, k++)
155  for (unsigned int j = 0; j < fe_collection[fe].dofs_per_cell; ++j)
156  fourier_transform_matrices[fe](k, j) =
157  integrate(fe_collection[fe],
158  q_collection[fe],
159  k_vectors(k1, k2, k3),
160  j);
161  }
162  }
163 } // namespace
164 
165 
166 
167 namespace FESeries
168 {
169  template <int dim, int spacedim>
171  const unsigned int N,
172  const hp::FECollection<dim, spacedim> &fe_collection,
173  const hp::QCollection<dim> & q_collection)
174  : fe_collection(&fe_collection)
175  , q_collection(&q_collection)
176  , fourier_transform_matrices(fe_collection.size())
177  {
178  set_k_vectors(k_vectors, N);
180  }
181 
182 
183 
184  template <int dim, int spacedim>
185  template <typename Number>
186  void
188  const Vector<Number> & local_dof_values,
189  const unsigned int cell_active_fe_index,
190  Table<dim, CoefficientType> &fourier_coefficients)
191  {
192  ensure_existence(*fe_collection,
193  *q_collection,
194  k_vectors,
195  cell_active_fe_index,
196  fourier_transform_matrices);
197 
198  const FullMatrix<CoefficientType> &matrix =
199  fourier_transform_matrices[cell_active_fe_index];
200 
201  std::fill(unrolled_coefficients.begin(),
202  unrolled_coefficients.end(),
203  CoefficientType(0.));
204 
205  Assert(unrolled_coefficients.size() == matrix.m(), ExcInternalError());
206 
207  Assert(local_dof_values.size() == matrix.n(),
208  ExcDimensionMismatch(local_dof_values.size(), matrix.n()));
209 
210  for (unsigned int i = 0; i < unrolled_coefficients.size(); i++)
211  for (unsigned int j = 0; j < local_dof_values.size(); j++)
212  unrolled_coefficients[i] += matrix[i][j] * local_dof_values[j];
213 
214  fourier_coefficients.fill(unrolled_coefficients.begin());
215  }
216 } // namespace FESeries
217 
218 
219 // explicit instantiations
220 #include "fe_series_fourier.inst"
221 
222 DEAL_II_NAMESPACE_CLOSE
void calculate(const ::Vector< Number > &local_dof_values, const unsigned int cell_active_fe_index, Table< dim, CoefficientType > &fourier_coefficients)
#define AssertIndexRange(index, range)
Definition: exceptions.h:1637
const Point< dim > & point(const unsigned int i) const
unsigned int size() const
size_type n_elements() const
#define Assert(cond, exc)
Definition: exceptions.h:1407
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
Fourier(const unsigned int size_in_each_direction, const hp::FECollection< dim, spacedim > &fe_collection, const hp::QCollection< dim > &q_collection)
SymmetricTensor< rank, dim, Number > sum(const SymmetricTensor< rank, dim, Number > &local, const MPI_Comm &mpi_communicator)
Table< dim, Tensor< 1, dim > > k_vectors
Definition: fe_series.h:132
unsigned int size() const
virtual double shape_value(const unsigned int i, const Point< dim > &p) const
Definition: fe.cc:159
Definition: mpi.h:90
static constexpr double PI
Definition: numbers.h:146
std::vector< CoefficientType > unrolled_coefficients
Definition: fe_series.h:142
Definition: table.h:37
void fill(InputIterator entries, const bool C_style_indexing=true)
double weight(const unsigned int i) const
static ::ExceptionBase & ExcInternalError()