90 const unsigned int face_no = 0;
94 std::vector<FullMatrix<double>> face_embeddings(
105 unsigned int target_row = 0;
106 for (
const auto &face_embedding : face_embeddings)
107 for (
unsigned int i = 0; i < face_embedding.m(); ++i)
109 for (
unsigned int j = 0; j < face_embedding.n(); ++j)
144 std::ostringstream namebuf;
146 namebuf <<
"FE_ABF<" << dim <<
">(" <<
rt_order <<
")";
148 return namebuf.str();
154std::unique_ptr<FiniteElement<dim, dim>>
157 return std::make_unique<FE_ABF<dim>>(
rt_order);
173 const unsigned int n_interior_points = cell_quadrature.size();
178 const unsigned int face_no = 0;
180 unsigned int n_face_points = (dim > 1) ? 1 : 0;
182 for (
unsigned int d = 1; d < dim; ++d)
183 n_face_points *= deg + 1;
192 std::array<std::unique_ptr<AnisotropicPolynomials<dim>>, dim> polynomials_abf;
195 for (
unsigned int dd = 0; dd < dim; ++dd)
197 std::vector<std::vector<Polynomials::Polynomial<double>>> poly(dim);
198 for (
unsigned int d = 0; d < dim; ++d)
202 polynomials_abf[dd] = std::make_unique<AnisotropicPolynomials<dim>>(poly);
206 unsigned int current = 0;
210 const QGauss<dim - 1> face_points(deg + 1);
219 for (
unsigned int k = 0; k < n_face_points; ++k)
222 face_points.point(k);
226 for (
unsigned int i = 0; i < legendre.n(); ++i)
229 face_points.weight(k) *
230 legendre.compute_value(i, face_points.point(k));
237 for (; current < GeometryInfo<dim>::faces_per_cell * n_face_points;
252 for (
unsigned int k = 0; k < faces.
size(); ++k)
254 for (
unsigned int i = 0; i < polynomials_abf[0]->n() * dim; ++i)
257 polynomials_abf[i % dim]->compute_value(i / dim,
268 std::array<std::unique_ptr<AnisotropicPolynomials<dim>>, dim> polynomials;
270 for (
unsigned int dd = 0; dd < dim; ++dd)
272 std::vector<std::vector<Polynomials::Polynomial<double>>> poly(dim);
273 for (
unsigned int d = 0; d < dim; ++d)
277 polynomials[dd] = std::make_unique<AnisotropicPolynomials<dim>>(poly);
283 for (
unsigned int k = 0; k < cell_quadrature.size(); ++k)
285 for (
unsigned int i = 0; i < polynomials[0]->n(); ++i)
286 for (
unsigned int d = 0; d < dim; ++d)
288 cell_quadrature.weight(k) *
289 polynomials[d]->compute_value(i, cell_quadrature.point(k));
296 for (
unsigned int k = 0; k < cell_quadrature.size(); ++k)
304 polynomials_abf[0]->n() * dim,
308 for (
unsigned int k = 0; k < cell_quadrature.size(); ++k)
310 for (
unsigned int i = 0; i < polynomials_abf[0]->n() * dim; ++i)
313 polynomials_abf[i % dim]->compute_grad(i / dim,
314 cell_quadrature.point(k)) *
315 cell_quadrature.weight(k);
318 for (
unsigned int d = 0; d < dim; ++d)
346 for (
unsigned int i = 0; i < GeometryInfo<dim>::max_children_per_cell;
353 const unsigned int n_face_points = q_base.size();
369 for (
unsigned int k = 0; k < q_face.
size(); ++k)
374 for (
unsigned int sub = 0; sub < GeometryInfo<dim>::max_children_per_face;
398 for (
unsigned int k = 0; k < n_face_points; ++k)
399 for (
unsigned int i_child = 0; i_child < this->n_dofs_per_cell();
401 for (
unsigned int i_face = 0;
411 face * this->n_dofs_per_face(face) + i_face, i_child) +=
413 cached_values_face(i_child, k) *
415 face * this->n_dofs_per_face(face) + i_face,
428 std::array<std::unique_ptr<AnisotropicPolynomials<dim>>, dim> polynomials;
429 for (
unsigned int dd = 0; dd < dim; ++dd)
431 std::vector<std::vector<Polynomials::Polynomial<double>>> poly(dim);
432 for (
unsigned int d = 0; d < dim; ++d)
436 polynomials[dd] = std::make_unique<AnisotropicPolynomials<dim>>(poly);
442 const unsigned int face_no = 0;
445 const unsigned int start_cell_dofs =
454 for (
unsigned int k = 0; k < q_cell.size(); ++k)
456 for (
unsigned int d = 0; d < dim; ++d)
457 cached_values_cell(i, k, d) =
460 for (
unsigned int child = 0; child < GeometryInfo<dim>::max_children_per_cell;
468 for (
unsigned int k = 0; k < q_sub.
size(); ++k)
469 for (
unsigned int i_child = 0; i_child < this->n_dofs_per_cell();
471 for (
unsigned int d = 0; d < dim; ++d)
472 for (
unsigned int i_weight = 0; i_weight < polynomials[d]->n();
475 this->
restriction[iso][child](start_cell_dofs + i_weight * dim +
478 q_sub.
weight(k) * cached_values_cell(i_child, k, d) *
479 polynomials[d]->compute_value(i_weight, q_sub.
point(k));
487std::vector<unsigned int>
493 return std::vector<unsigned int>();
502 for (
unsigned int d = 0; d < dim - 1; ++d)
508 std::vector<unsigned int> dpo(dim + 1);
510 dpo[dim] = interior_dofs;
522 const unsigned int face_index)
const
543 return (face_index !=
565 std::vector<double> &nodal_values)
const
567 Assert(support_point_values.size() == this->generalized_support_points.size(),
569 this->generalized_support_points.size()));
573 Assert(nodal_values.size() == this->n_dofs_per_cell(),
576 std::fill(nodal_values.begin(), nodal_values.end(), 0.);
580 for (
unsigned int k = 0; k < n_face_points; ++k)
585 support_point_values[face * n_face_points + k][
GeometryInfo<
592 const unsigned int face_no = 0;
594 const unsigned int start_cell_dofs =
596 const unsigned int start_cell_points =
601 for (
unsigned int d = 0; d < dim; ++d)
602 nodal_values[start_cell_dofs + i * dim + d] +=
604 support_point_values[k + start_cell_points][d];
606 const unsigned int start_abf_dofs =
612 for (
unsigned int d = 0; d < dim; ++d)
613 nodal_values[start_abf_dofs + i] +=
615 support_point_values[k + start_cell_points][d];
621 for (
unsigned int fp = 0; fp < n_face_points; ++fp)
631 nodal_values[start_abf_dofs + i] +=
633 support_point_values[face * n_face_points + fp][
GeometryInfo<
640 if (std::fabs(nodal_values[start_abf_dofs + i]) < 1.0e-16)
641 nodal_values[start_abf_dofs + i] = 0.0;
657#include "fe_abf.inst"
void initialize_quad_dof_index_permutation_and_sign_change()
Table< 2, double > boundary_weights
virtual std::unique_ptr< FiniteElement< dim, dim > > clone() 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
virtual std::size_t memory_consumption() const override
void initialize_restriction()
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 degree)
Table< 2, double > boundary_weights_abf
const unsigned int rt_order
void initialize_support_points(const unsigned int rt_degree)
Table< 3, double > interior_weights
Table< 3, double > interior_weights_abf
FullMatrix< double > inverse_node_matrix
virtual double shape_value_component(const unsigned int i, const Point< dim > &p, const unsigned int component) const override
std::vector< MappingKind > mapping_kind
FE_PolyTensor(const TensorPolynomialsBase< dim > &polynomials, const FiniteElementData< dim > &fe_data, const std::vector< bool > &restriction_is_additive_flags, const std::vector< ComponentMask > &nonzero_components)
unsigned int n_dofs_per_cell() const
const unsigned int dofs_per_face
unsigned int n_dofs_per_face(unsigned int face_no=0, unsigned int child=0) const
unsigned int n_components() const
ReferenceCell reference_cell() const
unsigned int n_unique_faces() const
FiniteElementData(const std::vector< unsigned int > &dofs_per_object, const unsigned int n_components, const unsigned int degree, const Conformity conformity=unknown, const BlockIndices &block_indices=BlockIndices())
std::vector< std::vector< FullMatrix< double > > > restriction
void reinit_restriction_and_prolongation_matrices(const bool isotropic_restriction_only=false, const bool isotropic_prolongation_only=false)
std::vector< std::vector< Point< dim - 1 > > > generalized_face_support_points
FullMatrix< double > interface_constraints
std::vector< Point< dim > > generalized_support_points
std::vector< std::vector< FullMatrix< double > > > prolongation
static std::vector< Polynomial< double > > generate_complete_basis(const unsigned int degree)
static std::vector< Polynomial< number > > generate_complete_basis(const unsigned int degree)
static DataSetDescriptor face(const ReferenceCell &reference_cell, const unsigned int face_no, const bool face_orientation, const bool face_flip, const bool face_rotation, const unsigned int n_quadrature_points)
static Quadrature< dim > project_to_child(const ReferenceCell &reference_cell, const Quadrature< dim > &quadrature, const unsigned int child_no)
static Quadrature< dim > project_to_all_faces(const ReferenceCell &reference_cell, const hp::QCollection< dim - 1 > &quadrature)
static void project_to_subface(const ReferenceCell &reference_cell, const SubQuadrature &quadrature, const unsigned int face_no, const unsigned int subface_no, std::vector< Point< dim > > &q_points, const RefinementCase< dim - 1 > &ref_case=RefinementCase< dim - 1 >::isotropic_refinement)
static void project_to_face(const ReferenceCell &reference_cell, const SubQuadrature &quadrature, const unsigned int face_no, std::vector< Point< dim > > &q_points)
const Point< dim > & point(const unsigned int i) const
double weight(const unsigned int i) const
unsigned int size() const
static constexpr unsigned char default_combined_face_orientation()
#define DEAL_II_NAMESPACE_OPEN
#define DEAL_II_NAMESPACE_CLOSE
constexpr std::array< unsigned int, GeometryInfo< dim >::faces_per_cell > GeometryInfo< dim >::unit_normal_direction
#define Assert(cond, exc)
static ::ExceptionBase & ExcImpossibleInDim(int arg1)
#define AssertDimension(dim1, dim2)
#define AssertIndexRange(index, range)
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
#define DEAL_II_NOT_IMPLEMENTED()
constexpr T fixed_power(const T t)
static unsigned int child_cell_on_face(const RefinementCase< dim > &ref_case, const unsigned int face, const unsigned int subface, const bool face_orientation=true, const bool face_flip=false, const bool face_rotation=false, const RefinementCase< dim - 1 > &face_refinement_case=RefinementCase< dim - 1 >::isotropic_refinement)
static constexpr std::array< int, faces_per_cell > unit_normal_orientation
static constexpr unsigned int faces_per_cell
static std_cxx20::ranges::iota_view< unsigned int, unsigned int > face_indices()
static constexpr std::array< unsigned int, faces_per_cell > unit_normal_direction
static constexpr std::array< unsigned int, faces_per_cell > opposite_face