33 for (
unsigned int i = 0; i <
N; ++i)
40 for (
unsigned int i = 0; i <
N; ++i)
41 for (
unsigned int j = 0; j <
N; ++j)
51 for (
unsigned int i = 0; i <
N; ++i)
52 for (
unsigned int j = 0; j <
N; ++j)
53 for (
unsigned int k = 0; k <
N; ++k)
63 template <
int dim,
int spacedim>
69 const unsigned int component)
71 std::complex<double>
sum = 0;
72 for (
unsigned int q = 0; q < quadrature.
size(); ++q)
75 sum +=
std::exp(std::complex<double>(0, 1) * (k_vector * x_q)) *
88 template <
int spacedim>
91 const std::vector<unsigned int> & n_coefficients_per_direction,
95 const unsigned int fe,
96 const unsigned int component,
97 std::vector<
FullMatrix<std::complex<double>>> &fourier_transform_matrices)
101 if (fourier_transform_matrices[fe].m() == 0)
103 fourier_transform_matrices[fe].reinit(
104 n_coefficients_per_direction[fe],
105 fe_collection[fe].n_dofs_per_cell());
107 for (
unsigned int k = 0; k < n_coefficients_per_direction[fe]; ++k)
108 for (
unsigned int j = 0; j < fe_collection[fe].n_dofs_per_cell(); ++j)
109 fourier_transform_matrices[fe](k, j) = integrate(
110 fe_collection[fe], q_collection[fe], k_vectors(k), j, component);
114 template <
int spacedim>
117 const std::vector<unsigned int> & n_coefficients_per_direction,
121 const unsigned int fe,
122 const unsigned int component,
123 std::vector<
FullMatrix<std::complex<double>>> &fourier_transform_matrices)
127 if (fourier_transform_matrices[fe].m() == 0)
129 fourier_transform_matrices[fe].reinit(
130 Utilities::fixed_power<2>(n_coefficients_per_direction[fe]),
131 fe_collection[fe].n_dofs_per_cell());
134 for (
unsigned int k1 = 0; k1 < n_coefficients_per_direction[fe]; ++k1)
135 for (
unsigned int k2 = 0; k2 < n_coefficients_per_direction[fe];
137 for (
unsigned int j = 0; j < fe_collection[fe].n_dofs_per_cell();
139 fourier_transform_matrices[fe](k, j) =
140 integrate(fe_collection[fe],
148 template <
int spacedim>
151 const std::vector<unsigned int> & n_coefficients_per_direction,
155 const unsigned int fe,
156 const unsigned int component,
157 std::vector<
FullMatrix<std::complex<double>>> &fourier_transform_matrices)
161 if (fourier_transform_matrices[fe].m() == 0)
163 fourier_transform_matrices[fe].reinit(
164 Utilities::fixed_power<3>(n_coefficients_per_direction[fe]),
165 fe_collection[fe].n_dofs_per_cell());
168 for (
unsigned int k1 = 0; k1 < n_coefficients_per_direction[fe]; ++k1)
169 for (
unsigned int k2 = 0; k2 < n_coefficients_per_direction[fe]; ++k2)
170 for (
unsigned int k3 = 0; k3 < n_coefficients_per_direction[fe];
172 for (
unsigned int j = 0; j < fe_collection[fe].n_dofs_per_cell();
174 fourier_transform_matrices[fe](k, j) =
175 integrate(fe_collection[fe],
177 k_vectors(k1, k2, k3),
188 template <
int dim,
int spacedim>
190 const std::vector<unsigned int> & n_coefficients_per_direction,
193 const unsigned int component_)
194 : n_coefficients_per_direction(n_coefficients_per_direction)
195 , fe_collection(&fe_collection)
196 , q_collection(q_collection)
197 , fourier_transform_matrices(fe_collection.size())
202 ExcMessage(
"All parameters are supposed to have the same size."));
208 "For vector-valued problems, you need to explicitly specify for "
209 "which vector component you will want to do a Fourier decomposition "
210 "by setting the 'component' argument of this constructor."));
214 const unsigned int max_n_coefficients_per_direction =
217 set_k_vectors(
k_vectors, max_n_coefficients_per_direction);
225 template <
int dim,
int spacedim>
227 const unsigned int n_coefficients_per_direction,
231 std::vector<unsigned
int>(fe_collection.size(),
232 n_coefficients_per_direction),
239 template <
int dim,
int spacedim>
255 template <
int dim,
int spacedim>
260 for (
unsigned int fe = 0; fe < fe_collection->
size(); ++fe)
262 ensure_existence(n_coefficients_per_direction,
268 fourier_transform_matrices);
276 template <
int dim,
int spacedim>
279 const unsigned int index)
const
281 return n_coefficients_per_direction[index];
286 template <
int dim,
int spacedim>
287 template <
typename Number>
291 const unsigned int cell_active_fe_index,
294 for (
unsigned int d = 0;
d < dim; ++
d)
296 n_coefficients_per_direction[cell_active_fe_index]);
298 ensure_existence(n_coefficients_per_direction,
302 cell_active_fe_index,
304 fourier_transform_matrices);
307 fourier_transform_matrices[cell_active_fe_index];
309 unrolled_coefficients.resize(Utilities::fixed_power<dim>(
310 n_coefficients_per_direction[cell_active_fe_index]));
311 std::fill(unrolled_coefficients.begin(),
312 unrolled_coefficients.end(),
320 for (
unsigned int i = 0; i < unrolled_coefficients.size(); i++)
321 for (
unsigned int j = 0; j < local_dof_values.
size(); j++)
322 unrolled_coefficients[i] +=
matrix[i][j] * local_dof_values[j];
324 fourier_coefficients.fill(unrolled_coefficients.begin());
330#include "fe_series_fourier.inst"
SmartPointer< const hp::FECollection< dim, spacedim > > fe_collection
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
std::vector< CoefficientType > unrolled_coefficients
const unsigned int component
typename std::complex< double > CoefficientType
unsigned int get_n_coefficients_per_direction(const unsigned int index) const
const std::vector< unsigned int > n_coefficients_per_direction
bool operator==(const Fourier< dim, spacedim > &fourier) const
std::vector< FullMatrix< CoefficientType > > fourier_transform_matrices
const hp::QCollection< dim > q_collection
Fourier(const std::vector< unsigned int > &n_coefficients_per_direction, const hp::FECollection< dim, spacedim > &fe_collection, const hp::QCollection< dim > &q_collection, const unsigned int component=numbers::invalid_unsigned_int)
void precalculate_all_transformation_matrices()
virtual double shape_value_component(const unsigned int i, const Point< dim > &p, const unsigned int component) const
const Point< dim > & point(const unsigned int i) const
double weight(const unsigned int i) const
unsigned int size() const
SymmetricTensor< rank, dim, Number > sum(const SymmetricTensor< rank, dim, Number > &local, const MPI_Comm &mpi_communicator)
unsigned int size() const
#define DEAL_II_NAMESPACE_OPEN
#define DEAL_II_NAMESPACE_CLOSE
#define Assert(cond, exc)
#define AssertDimension(dim1, dim2)
#define AssertIndexRange(index, range)
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
static ::ExceptionBase & ExcMessage(std::string arg1)
Task< RT > new_task(const std::function< RT()> &function)
@ matrix
Contents is actually a matrix.
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
static constexpr double PI
static const unsigned int invalid_unsigned_int
::VectorizedArray< Number, width > exp(const ::VectorizedArray< Number, width > &)