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;
185 this->generalized_support_points.resize(
187 this->generalized_face_support_points[face_no].resize(n_face_points);
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);
214 boundary_weights.reinit(n_face_points, legendre.n());
219 for (
unsigned int k = 0; k < n_face_points; ++k)
221 this->generalized_face_support_points[face_no][k] =
222 face_points.point(k);
226 for (
unsigned int i = 0; i < legendre.n(); ++i)
228 boundary_weights(k, i) =
229 face_points.weight(k) *
230 legendre.compute_value(i, face_points.point(k));
237 for (
unsigned int face_no = 0;
238 face_no < GeometryInfo<dim>::faces_per_cell;
242 this->reference_cell(),
246 for (
unsigned int face_point = 0; face_point < n_face_points;
250 this->generalized_support_points[current] =
251 faces.
point(offset + face_point);
262 boundary_weights_abf.reinit(faces.
size(), polynomials_abf[0]->n() * dim);
263 for (
unsigned int k = 0; k < faces.
size(); ++k)
265 for (
unsigned int i = 0; i < polynomials_abf[0]->n() * dim; ++i)
267 boundary_weights_abf(k, i) =
268 polynomials_abf[i % dim]->compute_value(i / dim,
279 std::array<std::unique_ptr<AnisotropicPolynomials<dim>>, dim> polynomials;
281 for (
unsigned int dd = 0; dd < dim; ++dd)
283 std::vector<std::vector<Polynomials::Polynomial<double>>> poly(dim);
284 for (
unsigned int d = 0; d < dim; ++d)
288 polynomials[dd] = std::make_unique<AnisotropicPolynomials<dim>>(poly);
291 interior_weights.reinit(
294 for (
unsigned int k = 0; k < cell_quadrature.size(); ++k)
296 for (
unsigned int i = 0; i < polynomials[0]->n(); ++i)
297 for (
unsigned int d = 0; d < dim; ++d)
298 interior_weights(k, i, d) =
299 cell_quadrature.weight(k) *
300 polynomials[d]->compute_value(i, cell_quadrature.point(k));
307 for (
unsigned int k = 0; k < cell_quadrature.size(); ++k)
308 this->generalized_support_points[current++] = cell_quadrature.point(k);
315 polynomials_abf[0]->n() * dim,
319 for (
unsigned int k = 0; k < cell_quadrature.size(); ++k)
321 for (
unsigned int i = 0; i < polynomials_abf[0]->n() * dim; ++i)
324 polynomials_abf[i % dim]->compute_grad(i / dim,
325 cell_quadrature.point(k)) *
326 cell_quadrature.weight(k);
329 for (
unsigned int d = 0; d < dim; ++d)
330 interior_weights_abf(k, i, d) = -poly_grad[d];
334 Assert(current == this->generalized_support_points.size(),
357 for (
unsigned int i = 0; i < GeometryInfo<dim>::max_children_per_cell;
359 this->restriction[iso][i].reinit(0, 0);
363 const QGauss<dim - 1> q_base(rt_order + 1);
364 const unsigned int n_face_points = q_base.size();
374 this->reference_cell(),
383 for (
unsigned int k = 0; k < q_face.
size(); ++k)
384 for (
unsigned int i = 0; i < this->n_dofs_per_cell(); ++i)
385 cached_values_face(i, k) = this->shape_value_component(
388 for (
unsigned int sub = 0; sub < GeometryInfo<dim>::max_children_per_face;
396 this->reference_cell(),
417 for (
unsigned int k = 0; k < n_face_points; ++k)
418 for (
unsigned int i_child = 0; i_child < this->n_dofs_per_cell();
420 for (
unsigned int i_face = 0;
421 i_face < this->n_dofs_per_face(face);
429 this->restriction[iso][child](
430 face * this->n_dofs_per_face(face) + i_face, i_child) +=
431 Utilities::fixed_power<dim - 1>(.5) * q_sub.
weight(k) *
432 cached_values_face(i_child, k) *
433 this->shape_value_component(
434 face * this->n_dofs_per_face(face) + i_face,
447 std::array<std::unique_ptr<AnisotropicPolynomials<dim>>, dim> polynomials;
448 for (
unsigned int dd = 0; dd < dim; ++dd)
450 std::vector<std::vector<Polynomials::Polynomial<double>>> poly(dim);
451 for (
unsigned int d = 0; d < dim; ++d)
455 polynomials[dd] = std::make_unique<AnisotropicPolynomials<dim>>(poly);
461 const unsigned int face_no = 0;
464 const unsigned int start_cell_dofs =
473 for (
unsigned int k = 0; k < q_cell.size(); ++k)
474 for (
unsigned int i = 0; i < this->n_dofs_per_cell(); ++i)
475 for (
unsigned int d = 0; d < dim; ++d)
476 cached_values_cell(i, k, d) =
477 this->shape_value_component(i, q_cell.point(k), d);
479 for (
unsigned int child = 0; child < GeometryInfo<dim>::max_children_per_cell;
487 for (
unsigned int k = 0; k < q_sub.
size(); ++k)
488 for (
unsigned int i_child = 0; i_child < this->n_dofs_per_cell();
490 for (
unsigned int d = 0; d < dim; ++d)
491 for (
unsigned int i_weight = 0; i_weight < polynomials[d]->n();
494 this->restriction[iso][child](start_cell_dofs + i_weight * dim +
497 q_sub.
weight(k) * cached_values_cell(i_child, k, d) *
498 polynomials[d]->compute_value(i_weight, q_sub.
point(k));
506std::vector<unsigned int>
512 return std::vector<unsigned int>();
520 unsigned int dofs_per_face = 1;
521 for (
unsigned int d = 0; d < dim - 1; ++d)
522 dofs_per_face *= rt_order + 1;
525 const unsigned int interior_dofs = dim * (rt_order + 1) * dofs_per_face;
527 std::vector<unsigned int> dpo(dim + 1);
528 dpo[dim - 1] = dofs_per_face;
529 dpo[dim] = interior_dofs;
541 const unsigned int face_index)
const
562 return (face_index !=
584 std::vector<double> &nodal_values)
const
586 Assert(support_point_values.size() == this->generalized_support_points.size(),
588 this->generalized_support_points.size()));
589 Assert(support_point_values[0].
size() == this->n_components(),
591 this->n_components()));
592 Assert(nodal_values.size() == this->n_dofs_per_cell(),
595 std::fill(nodal_values.begin(), nodal_values.end(), 0.);
597 const unsigned int n_face_points = boundary_weights.size(0);
599 for (
unsigned int k = 0; k < n_face_points; ++k)
600 for (
unsigned int i = 0; i < boundary_weights.size(1); ++i)
602 nodal_values[i + face * this->n_dofs_per_face(face)] +=
603 boundary_weights(k, i) *
604 support_point_values[face * n_face_points + k][
GeometryInfo<
605 dim>::unit_normal_direction[face]];
611 const unsigned int face_no = 0;
613 const unsigned int start_cell_dofs =
615 const unsigned int start_cell_points =
618 for (
unsigned int k = 0; k < interior_weights.size(0); ++k)
619 for (
unsigned int i = 0; i < interior_weights.size(1); ++i)
620 for (
unsigned int d = 0; d < dim; ++d)
621 nodal_values[start_cell_dofs + i * dim + d] +=
622 interior_weights(k, i, d) *
623 support_point_values[k + start_cell_points][d];
625 const unsigned int start_abf_dofs =
626 start_cell_dofs + interior_weights.size(1) * dim;
629 for (
unsigned int k = 0; k < interior_weights_abf.size(0); ++k)
630 for (
unsigned int i = 0; i < interior_weights_abf.size(1); ++i)
631 for (
unsigned int d = 0; d < dim; ++d)
632 nodal_values[start_abf_dofs + i] +=
633 interior_weights_abf(k, i, d) *
634 support_point_values[k + start_cell_points][d];
640 for (
unsigned int fp = 0; fp < n_face_points; ++fp)
645 this->reference_cell(),
649 for (
unsigned int i = 0; i < boundary_weights_abf.size(1); ++i)
650 nodal_values[start_abf_dofs + i] +=
651 n_orient * boundary_weights_abf(k + fp, i) *
652 support_point_values[face * n_face_points + fp][
GeometryInfo<
653 dim>::unit_normal_direction[face]];
658 for (
unsigned int i = 0; i < boundary_weights_abf.size(1); ++i)
659 if (std::fabs(nodal_values[start_abf_dofs + i]) < 1.0e-16)
660 nodal_values[start_abf_dofs + i] = 0.0;
676#include "fe/fe_abf.inst"
void initialize_quad_dof_index_permutation_and_sign_change()
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)
void initialize_support_points(const unsigned int rt_degree)
FullMatrix< double > inverse_node_matrix
std::vector< MappingKind > mapping_kind
unsigned int n_dofs_per_cell() const
unsigned int n_dofs_per_face(unsigned int face_no=0, unsigned int child=0) const
unsigned int n_unique_faces() const
void reinit_restriction_and_prolongation_matrices(const bool isotropic_restriction_only=false, const bool isotropic_prolongation_only=false)
FullMatrix< double > interface_constraints
std::vector< std::vector< FullMatrix< double > > > prolongation
void invert(const FullMatrix< number2 > &M)
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_all_faces(const ReferenceCell &reference_cell, const hp::QCollection< dim - 1 > &quadrature)
static Quadrature< dim > project_to_child(const ReferenceCell &reference_cell, const Quadrature< dim > &quadrature, const unsigned int child_no)
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
#define DEAL_II_NAMESPACE_OPEN
#define DEAL_II_NAMESPACE_CLOSE
#define DEAL_II_NOT_IMPLEMENTED()
#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)
constexpr types::geometric_orientation default_geometric_orientation
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 std_cxx20::ranges::iota_view< unsigned int, unsigned int > face_indices()