17 #include <deal.II/base/quadrature_lib.h> 18 #include <deal.II/base/std_cxx14/memory.h> 20 #include <deal.II/fe/fe_face.h> 21 #include <deal.II/fe/fe_nothing.h> 22 #include <deal.II/fe/fe_poly_face.templates.h> 23 #include <deal.II/fe/fe_tools.h> 25 #include <deal.II/lac/householder.h> 29 DEAL_II_NAMESPACE_OPEN
34 namespace FE_FaceQImplementation
39 get_QGaussLobatto_points(
const unsigned int degree)
44 return std::vector<Point<1>>(1,
Point<1>(0.5));
50 template <
int dim,
int spacedim>
55 internal::FE_FaceQImplementation::get_QGaussLobatto_points(degree))),
60 std::vector<bool>(1, true))
63 const unsigned int codim = dim - 1;
65 Utilities::fixed_power<codim>(this->degree + 1));
67 if (this->degree == 0)
68 for (
unsigned int d = 0; d < codim; ++d)
72 std::vector<Point<1>> points =
73 internal::FE_FaceQImplementation::get_QGaussLobatto_points(degree);
76 for (
unsigned int iz = 0; iz <= ((codim > 2) ? this->degree : 0); ++iz)
77 for (
unsigned int iy = 0; iy <= ((codim > 1) ? this->degree : 0); ++iy)
78 for (
unsigned int ix = 0; ix <= this->
degree; ++ix)
98 for (
unsigned int i = 0; i < n_face_dofs; ++i)
99 for (
unsigned int d = 0; d < dim; ++d)
101 for (
unsigned int e = 0, c = 0; e < dim; ++e)
105 unsigned int renumber = i;
106 if (dim == 3 && d == 1)
107 renumber = i / (degree + 1) + (degree + 1) * (i % (degree + 1));
120 template <
int dim,
int spacedim>
121 std::unique_ptr<FiniteElement<dim, spacedim>>
124 return std_cxx14::make_unique<FE_FaceQ<dim, spacedim>>(this->
degree);
129 template <
int dim,
int spacedim>
136 std::ostringstream namebuf;
140 return namebuf.str();
145 template <
int dim,
int spacedim>
153 interpolation_matrix);
158 template <
int dim,
int spacedim>
162 const unsigned int subface,
189 source_fe->get_unit_face_support_points());
194 const double eps = 2e-13 * (this->
degree + 1) * (dim - 1);
198 for (
unsigned int i = 0; i < source_fe->dofs_per_face; ++i)
200 const Point<dim - 1> p =
202 face_quadrature.point(i) :
204 face_quadrature.point(i), subface);
213 if (std::fabs(matrix_entry - 1.0) < eps)
215 if (std::fabs(matrix_entry) < eps)
218 interpolation_matrix(i, j) = matrix_entry;
224 for (
unsigned int j = 0; j < source_fe->dofs_per_face; ++j)
229 sum += interpolation_matrix(j, i);
234 else if (
dynamic_cast<const FE_Nothing<dim> *
>(&x_source_fe) !=
nullptr)
247 template <
int dim,
int spacedim>
250 const unsigned int shape_index,
251 const unsigned int face_index)
const 258 template <
int dim,
int spacedim>
259 std::vector<unsigned int>
262 std::vector<unsigned int> dpo(dim + 1, 0U);
263 dpo[dim - 1] = deg + 1;
264 for (
unsigned int i = 1; i < dim - 1; ++i)
265 dpo[dim - 1] *= deg + 1;
271 template <
int dim,
int spacedim>
280 template <
int dim,
int spacedim>
281 std::vector<std::pair<unsigned int, unsigned int>>
286 return std::vector<std::pair<unsigned int, unsigned int>>();
291 template <
int dim,
int spacedim>
292 std::vector<std::pair<unsigned int, unsigned int>>
300 return std::vector<std::pair<unsigned int, unsigned int>>();
313 const unsigned int p = this->
degree;
314 const unsigned int q = fe_q_other->degree;
316 std::vector<std::pair<unsigned int, unsigned int>> identities;
318 const std::vector<unsigned int> &index_map_inverse =
320 const std::vector<unsigned int> &index_map_inverse_other =
321 fe_q_other->poly_space.get_numbering_inverse();
323 for (
unsigned int i = 0; i < p + 1; ++i)
324 for (
unsigned int j = 0; j < q + 1; ++j)
327 fe_q_other->unit_support_points[index_map_inverse_other[j]]
329 identities.emplace_back(i, j);
337 return std::vector<std::pair<unsigned int, unsigned int>>();
348 return std::vector<std::pair<unsigned int, unsigned int>>();
353 return std::vector<std::pair<unsigned int, unsigned int>>();
360 template <
int dim,
int spacedim>
361 std::vector<std::pair<unsigned int, unsigned int>>
369 return std::vector<std::pair<unsigned int, unsigned int>>();
381 const unsigned int p = this->
degree;
382 const unsigned int q = fe_q_other->degree;
384 std::vector<std::pair<unsigned int, unsigned int>> identities;
386 const std::vector<unsigned int> &index_map_inverse =
388 const std::vector<unsigned int> &index_map_inverse_other =
389 fe_q_other->poly_space.get_numbering_inverse();
391 std::vector<std::pair<unsigned int, unsigned int>> identities_1d;
393 for (
unsigned int i = 0; i < p + 1; ++i)
394 for (
unsigned int j = 0; j < q + 1; ++j)
397 fe_q_other->unit_support_points[index_map_inverse_other[j]]
399 identities_1d.emplace_back(i, j);
401 for (
unsigned int n1 = 0; n1 < identities_1d.size(); ++n1)
402 for (
unsigned int n2 = 0; n2 < identities_1d.size(); ++n2)
403 identities.emplace_back(identities_1d[n1].first * (p + 1) +
404 identities_1d[n2].first,
405 identities_1d[n1].second * (q + 1) +
406 identities_1d[n2].second);
414 return std::vector<std::pair<unsigned int, unsigned int>>();
425 return std::vector<std::pair<unsigned int, unsigned int>>();
430 return std::vector<std::pair<unsigned int, unsigned int>>();
437 template <
int dim,
int spacedim>
445 if (this->degree < fe_q_other->
degree)
447 else if (this->degree == fe_q_other->degree)
455 if (fe_nothing->is_dominating())
474 template <
int dim,
int spacedim>
475 std::pair<Table<2, bool>, std::vector<unsigned int>>
480 constant_modes(0, i) =
true;
481 return std::pair<Table<2, bool>, std::vector<unsigned int>>(
482 constant_modes, std::vector<unsigned int>(1, 0));
485 template <
int dim,
int spacedim>
489 std::vector<double> & nodal_values)
const 500 nodal_values[i] = support_point_values[i](0);
506 template <
int spacedim>
513 std::vector<bool>(1, true),
526 template <
int spacedim>
527 std::unique_ptr<FiniteElement<1, spacedim>>
530 return std_cxx14::make_unique<FE_FaceQ<1, spacedim>>(this->
degree);
535 template <
int spacedim>
542 std::ostringstream namebuf;
546 return namebuf.str();
551 template <
int spacedim>
559 interpolation_matrix);
564 template <
int spacedim>
577 interpolation_matrix(0, 0) = 1.;
582 template <
int spacedim>
585 const unsigned int face_index)
const 588 return (face_index == shape_index);
593 template <
int spacedim>
594 std::vector<unsigned int>
597 std::vector<unsigned int> dpo(2, 0U);
604 template <
int spacedim>
611 template <
int spacedim>
612 std::vector<std::pair<unsigned int, unsigned int>>
617 return std::vector<std::pair<unsigned int, unsigned int>>(1,
624 template <
int spacedim>
625 std::vector<std::pair<unsigned int, unsigned int>>
630 return std::vector<std::pair<unsigned int, unsigned int>>();
635 template <
int spacedim>
636 std::vector<std::pair<unsigned int, unsigned int>>
641 return std::vector<std::pair<unsigned int, unsigned int>>();
646 template <
int spacedim>
656 template <
int spacedim>
657 std::pair<Table<2, bool>, std::vector<unsigned int>>
662 constant_modes(0, i) =
true;
663 return std::pair<Table<2, bool>, std::vector<unsigned int>>(
664 constant_modes, std::vector<unsigned int>(1, 0));
669 template <
int spacedim>
685 template <
int spacedim>
693 const ::internal::FEValuesImplementation::MappingRelatedData<1,
706 template <
int spacedim>
710 const unsigned int face,
714 const ::internal::FEValuesImplementation::MappingRelatedData<1,
722 const unsigned int foffset = face;
726 output_data.shape_values(k, 0) = 0.;
727 output_data.shape_values(foffset, 0) = 1;
732 template <
int spacedim>
741 const ::internal::FEValuesImplementation::MappingRelatedData<1,
756 template <
int dim,
int spacedim>
760 Polynomials::Legendre::generate_complete_basis(degree)),
765 std::vector<bool>(1, true))
770 template <
int dim,
int spacedim>
771 std::unique_ptr<FiniteElement<dim, spacedim>>
774 return std_cxx14::make_unique<FE_FaceP<dim, spacedim>>(this->
degree);
779 template <
int dim,
int spacedim>
786 std::ostringstream namebuf;
790 return namebuf.str();
795 template <
int dim,
int spacedim>
798 const unsigned int shape_index,
799 const unsigned int face_index)
const 806 template <
int dim,
int spacedim>
807 std::vector<unsigned int>
810 std::vector<unsigned int> dpo(dim + 1, 0U);
811 dpo[dim - 1] = deg + 1;
812 for (
unsigned int i = 1; i < dim - 1; ++i)
814 dpo[dim - 1] *= deg + 1 + i;
815 dpo[dim - 1] /= i + 1;
822 template <
int dim,
int spacedim>
831 template <
int dim,
int spacedim>
839 if (this->degree < fe_q_other->
degree)
841 else if (this->degree == fe_q_other->degree)
849 if (fe_nothing->is_dominating())
868 template <
int dim,
int spacedim>
876 interpolation_matrix);
881 template <
int dim,
int spacedim>
885 const unsigned int subface,
913 const QGauss<dim - 1> face_quadrature(source_fe->degree + 1);
918 const double eps = 2e-13 * (this->
degree + 1) * (dim - 1);
922 for (
unsigned int k = 0; k < face_quadrature.size(); ++k)
924 const Point<dim - 1> p =
926 face_quadrature.point(k) :
928 face_quadrature.point(k), subface);
930 for (
unsigned int j = 0; j < source_fe->dofs_per_face; ++j)
931 mass(k, j) = source_fe->poly_space.compute_value(j, p);
943 for (
unsigned int k = 0; k < face_quadrature.size(); ++k)
945 const Point<dim - 1> p =
947 face_quadrature.point(k) :
949 face_quadrature.point(k), subface);
956 for (
unsigned int j = 0; j < source_fe->dofs_per_face; ++j)
958 double matrix_entry = v_out(j);
963 if (std::fabs(matrix_entry - 1.0) < eps)
965 if (std::fabs(matrix_entry) < eps)
968 interpolation_matrix(j, i) = matrix_entry;
972 else if (
dynamic_cast<const FE_Nothing<dim> *
>(&x_source_fe) !=
nullptr)
985 template <
int dim,
int spacedim>
986 std::pair<Table<2, bool>, std::vector<unsigned int>>
990 for (
unsigned int face = 0; face < GeometryInfo<dim>::faces_per_cell; ++face)
992 return std::pair<Table<2, bool>, std::vector<unsigned int>>(
993 constant_modes, std::vector<unsigned int>(1, 0));
998 template <
int spacedim>
1005 template <
int spacedim>
1012 std::ostringstream namebuf;
1016 return namebuf.str();
1022 #include "fe_face.inst" 1025 DEAL_II_NAMESPACE_CLOSE
static Point< dim > child_to_cell_coordinates(const Point< dim > &p, const unsigned int child_index, const RefinementCase< dim > refine_case=RefinementCase< dim >::isotropic_refinement)
Transformed quadrature weights.
virtual void get_face_interpolation_matrix(const FiniteElement< dim, spacedim > &source, FullMatrix< double > &matrix) const override
static const unsigned int invalid_unsigned_int
static::ExceptionBase & ExcLeastSquaresError(double arg1)
virtual UpdateFlags requires_update_flags(const UpdateFlags update_flags) const override
virtual std::pair< Table< 2, bool >, std::vector< unsigned int > > get_constant_modes() const override
const std::vector< Point< dim > > & get_unit_support_points() const
#define AssertDimension(dim1, dim2)
double compute_value(const unsigned int i, const Point< dim > &p) const
#define AssertIndexRange(index, range)
const unsigned int degree
std::vector< Point< dim-1 > > unit_face_support_points
virtual bool has_support_on_face(const unsigned int shape_index, const unsigned int face_index) const override
#define AssertThrow(cond, exc)
virtual std::unique_ptr< FiniteElement< dim, spacedim > > clone() const override
static::ExceptionBase & ExcInterpolationNotImplemented()
virtual std::unique_ptr< FiniteElement< dim, spacedim > > clone() const override
const std::vector< unsigned int > & get_numbering_inverse() const
std::vector< Point< dim > > unit_support_points
static::ExceptionBase & ExcMessage(std::string arg1)
virtual void get_subface_interpolation_matrix(const FiniteElement< dim, spacedim > &source, const unsigned int subface, FullMatrix< double > &matrix) const override
virtual std::vector< std::pair< unsigned int, unsigned int > > hp_line_dof_identities(const FiniteElement< dim, spacedim > &fe_other) const override
#define Assert(cond, exc)
static::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
Abstract base class for mapping classes.
virtual void get_subface_interpolation_matrix(const FiniteElement< dim, spacedim > &source, const unsigned int subface, FullMatrix< double > &matrix) const override
virtual FiniteElementDomination::Domination compare_for_face_domination(const FiniteElement< dim, spacedim > &fe_other) const override
FE_FaceQ(const unsigned int p)
virtual std::vector< std::pair< unsigned int, unsigned int > > hp_vertex_dof_identities(const FiniteElement< dim, spacedim > &fe_other) const override
virtual bool has_support_on_face(const unsigned int shape_index, const unsigned int face_index) const override
virtual std::string get_name() const override
static std::vector< unsigned int > get_dpo_vector(const unsigned int deg)
virtual bool hp_constraints_are_implemented() const override
Second derivatives of shape functions.
std::string dim_string(const int dim, const int spacedim)
virtual std::vector< std::pair< unsigned int, unsigned int > > hp_quad_dof_identities(const FiniteElement< dim, spacedim > &fe_other) const override
TensorProductPolynomials< dim-1 > poly_space
const unsigned int dofs_per_cell
Shape function gradients.
static std::vector< unsigned int > get_dpo_vector(const unsigned int deg)
const unsigned int dofs_per_face
virtual std::pair< Table< 2, bool >, std::vector< unsigned int > > get_constant_modes() const override
static::ExceptionBase & ExcNotImplemented()
virtual void get_face_interpolation_matrix(const FiniteElement< dim, spacedim > &source, FullMatrix< double > &matrix) const override
virtual FiniteElementDomination::Domination compare_for_face_domination(const FiniteElement< dim, spacedim > &fe_other) const override
double least_squares(Vector< number2 > &dst, const Vector< number2 > &src) const
virtual std::string get_name() const override
virtual void convert_generalized_support_point_values_to_dof_values(const std::vector< Vector< double >> &support_point_values, std::vector< double > &nodal_values) const override
double compute_value(const unsigned int i, const Point< dim > &p) const
Covariant transformation.
static::ExceptionBase & ExcInternalError()
virtual bool hp_constraints_are_implemented() const override