18 #include <deal.II/fe/fe_series.h> 19 #include <deal.II/base/numbers.h> 24 #include <deal.II/base/config.h> 25 #ifdef DEAL_II_WITH_GSL 26 #include <gsl/gsl_sf_legendre.h> 30 DEAL_II_NAMESPACE_OPEN
41 for (
unsigned int i=0; i<N; ++i)
50 for (
unsigned int i=0; i<N; ++i)
51 for (
unsigned int j=0; j<N; ++j)
62 for (
unsigned int i=0; i<N; ++i)
63 for (
unsigned int j=0; j<N; ++j)
64 for (
unsigned int k=0; k<N; ++k)
78 fe_collection(&fe_collection),
79 q_collection(&q_collection),
80 fourier_transform_matrices(fe_collection.size())
88 const unsigned int cell_active_fe_index,
89 Table<dim,std::complex<double> > &fourier_coefficients)
96 std::complex<double>(0.));
105 for (
unsigned int j = 0; j < local_dof_values.
size(); j++)
116 const unsigned int j)
118 std::complex<double> sum = 0;
119 for (
unsigned int q=0; q<quadrature.
size(); ++q)
122 sum += std::exp(std::complex<double>(0,1) *
134 Assert (fe < fe_collection->size(),
140 (*fe_collection)[fe].dofs_per_cell);
145 for (
unsigned int j=0; j<(*fe_collection)[fe].dofs_per_cell; ++j)
156 Assert (fe < fe_collection->size(),
162 (*fe_collection)[fe].dofs_per_cell);
168 for (
unsigned int j=0; j<(*fe_collection)[fe].dofs_per_cell; ++j)
179 Assert (fe < fe_collection->size(),
185 (*fe_collection)[fe].dofs_per_cell);
188 for (
unsigned int j=0; j<(*fe_collection)[fe].dofs_per_cell; ++j)
199 <<
"x["<<arg1 <<
"] = "<<arg2 <<
" is not in [0,1]");
208 #ifdef DEAL_II_WITH_GSL 210 for (
unsigned int d = 0; d < dim; d++)
212 const double x = 2.0*(x_q[d]-0.5);
213 Assert ( (x_q[d] <= 1.0) && (x_q[d] >= 0.),
215 const int ind = indices[d];
216 res *= sqrt(2.0) * gsl_sf_legendre_Pl (ind, x);
225 "in order to use Legendre transformation."));
237 for (
unsigned int d = 0; d < dim; d++)
238 res *= (0.5+indices[d]);
249 N(size_in_each_direction),
250 fe_collection(&fe_collection),
251 q_collection(&q_collection),
252 legendre_transform_matrices(fe_collection.size()),
260 const unsigned int cell_active_fe_index,
273 Assert (local_dof_values.size() == matrix.
n(),
277 for (
unsigned int j = 0; j < local_dof_values.size(); j++)
289 const unsigned int dof)
292 for (
unsigned int q=0; q<quadrature.
size(); ++q)
295 sum += Lh(x_q, indices) *
299 return sum * multiplier(indices);
305 Assert (fe < fe_collection->size(),
312 for (
unsigned int k=0; k<
N; ++k)
313 for (
unsigned int j=0; j<(*fe_collection)[fe].dofs_per_cell; ++j)
327 Assert (fe < fe_collection->size(),
336 for (
unsigned int k1=0; k1<
N; ++k1)
337 for (
unsigned int k2=0; k2<
N; ++k2,k++)
338 for (
unsigned int j=0; j<(*fe_collection)[fe].dofs_per_cell; ++j)
350 Assert (fe < fe_collection->size(),
359 for (
unsigned int k1=0; k1<
N; ++k1)
360 for (
unsigned int k2=0; k2<
N; ++k2)
361 for (
unsigned int k3=0; k3<
N; ++k3, k++)
362 for (
unsigned int j=0; j<(*fe_collection)[fe].dofs_per_cell; ++j)
373 const std::vector<double> &y)
378 Assert (x.size() == y.size(),
379 ExcMessage(
"x and y are expected to have the same size"));
382 ::ExcMessage(
"at least two points are required for linear regression fit"));
390 for (
unsigned int i = 0; i < x.size(); i++)
408 invK.
vmult(X,B,
false);
410 return std::make_pair(X(1),X(0));
428 DEAL_II_NAMESPACE_CLOSE
virtual double shape_value(const unsigned int i, const Point< dim > &p) const
static::ExceptionBase & ExcLegendre(int arg1, double arg2)
#define DeclException2(Exception2, type1, type2, outsequence)
void ensure_existence(const unsigned int fe_index)
SmartPointer< const hp::FECollection< dim > > fe_collection
Legendre(const unsigned int size_in_each_direction, const hp::FECollection< dim > &fe_collection, const hp::QCollection< dim > &q_collection)
std::pair< double, double > linear_regression(const std::vector< double > &x, const std::vector< double > &y)
SmartPointer< const hp::QCollection< dim > > q_collection
unsigned int size() const
void vmult(Vector< number2 > &w, const Vector< number2 > &v, const bool adding=false) const
#define AssertThrow(cond, exc)
void calculate(const ::Vector< double > &local_dof_values, const unsigned int cell_active_fe_index, Table< dim, double > &legendre_coefficients)
Fourier(const unsigned int size_in_each_direction, const hp::FECollection< dim > &fe_collection, const hp::QCollection< dim > &q_collection)
static::ExceptionBase & ExcIndexRange(int arg1, int arg2, int arg3)
const Point< dim > & point(const unsigned int i) const
std::vector< std::complex< double > > unrolled_coefficients
SmartPointer< const hp::QCollection< dim > > q_collection
unsigned int size() const
static::ExceptionBase & ExcMessage(std::string arg1)
void invert(const FullMatrix< number2 > &M)
#define Assert(cond, exc)
static::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
double weight(const unsigned int i) const
SmartPointer< const hp::FECollection< dim > > fe_collection
void calculate(const ::Vector< double > &local_dof_values, const unsigned int cell_active_fe_index, Table< dim, std::complex< double > > &fourier_coefficients)
void ensure_existence(const unsigned int fe_index)
size_type n_elements() const
std::vector< double > unrolled_coefficients
Table< dim, Tensor< 1, dim > > k_vectors
unsigned int size(const unsigned int i) const
std::vector< FullMatrix< std::complex< double > > > fourier_transform_matrices
void fill(InputIterator entries, const bool C_style_indexing=true)
std::vector< FullMatrix< double > > legendre_transform_matrices
static::ExceptionBase & ExcInternalError()