17 #include <deal.II/base/qprojector.h> 18 #include <deal.II/base/quadrature.h> 19 #include <deal.II/base/quadrature_lib.h> 20 #include <deal.II/base/std_cxx14/memory.h> 21 #include <deal.II/base/table.h> 22 #include <deal.II/base/utilities.h> 24 #include <deal.II/dofs/dof_accessor.h> 26 #include <deal.II/fe/fe.h> 27 #include <deal.II/fe/fe_abf.h> 28 #include <deal.II/fe/fe_tools.h> 29 #include <deal.II/fe/fe_values.h> 30 #include <deal.II/fe/mapping.h> 32 #include <deal.II/grid/tria.h> 33 #include <deal.II/grid/tria_iterator.h> 44 DEAL_II_NAMESPACE_OPEN
57 std::vector<bool>(dim, true)))
93 std::vector<FullMatrix<double>> face_embeddings(
102 unsigned int target_row = 0;
103 for (
unsigned int d = 0; d < face_embeddings.size(); ++d)
104 for (
unsigned int i = 0; i < face_embeddings[d].m(); ++i)
106 for (
unsigned int j = 0; j < face_embeddings[d].n(); ++j)
125 std::ostringstream namebuf;
127 namebuf <<
"FE_ABF<" << dim <<
">(" <<
rt_order <<
")";
129 return namebuf.str();
135 std::unique_ptr<FiniteElement<dim, dim>>
138 return std_cxx14::make_unique<FE_ABF<dim>>(
rt_order);
154 const unsigned int n_interior_points = cell_quadrature.
size();
156 unsigned int n_face_points = (dim > 1) ? 1 : 0;
158 for (
unsigned int d = 1; d < dim; ++d)
159 n_face_points *= deg + 1;
168 std::vector<AnisotropicPolynomials<dim> *> polynomials_abf(dim);
171 for (
unsigned int dd = 0; dd < dim; ++dd)
173 std::vector<std::vector<Polynomials::Polynomial<double>>> poly(dim);
174 for (
unsigned int d = 0; d < dim; ++d)
182 unsigned int current = 0;
186 QGauss<dim - 1> face_points(deg + 1);
195 for (
unsigned int k = 0; k < n_face_points; ++k)
201 for (
unsigned int i = 0; i < legendre.n(); ++i)
204 face_points.weight(k) *
205 legendre.compute_value(i, face_points.point(k));
211 for (; current < GeometryInfo<dim>::faces_per_cell * n_face_points;
226 for (
unsigned int k = 0; k < faces.
size(); ++k)
228 for (
unsigned int i = 0; i < polynomials_abf[0]->n() * dim; ++i)
231 polynomials_abf[i % dim]->compute_value(i / dim,
242 std::vector<AnisotropicPolynomials<dim> *> polynomials(dim);
244 for (
unsigned int dd = 0; dd < dim; ++dd)
246 std::vector<std::vector<Polynomials::Polynomial<double>>> poly(dim);
247 for (
unsigned int d = 0; d < dim; ++d)
257 for (
unsigned int k = 0; k < cell_quadrature.
size(); ++k)
259 for (
unsigned int i = 0; i < polynomials[0]->n(); ++i)
260 for (
unsigned int d = 0; d < dim; ++d)
262 cell_quadrature.
weight(k) *
263 polynomials[d]->compute_value(i, cell_quadrature.
point(k));
266 for (
unsigned int d = 0; d < dim; ++d)
267 delete polynomials[d];
273 for (
unsigned int k = 0; k < cell_quadrature.
size(); ++k)
281 polynomials_abf[0]->n() * dim,
285 for (
unsigned int k = 0; k < cell_quadrature.
size(); ++k)
287 for (
unsigned int i = 0; i < polynomials_abf[0]->n() * dim; ++i)
290 polynomials_abf[i % dim]->compute_grad(i / dim,
291 cell_quadrature.
point(k)) *
292 cell_quadrature.
weight(k);
295 for (
unsigned int d = 0; d < dim; ++d)
300 for (
unsigned int d = 0; d < dim; ++d)
301 delete polynomials_abf[d];
326 for (
unsigned int i = 0; i < GeometryInfo<dim>::max_children_per_cell;
333 const unsigned int n_face_points = q_base.size();
336 for (
unsigned int face = 0; face < GeometryInfo<dim>::faces_per_cell; ++face)
347 for (
unsigned int k = 0; k < q_face.
size(); ++k)
352 for (
unsigned int sub = 0; sub < GeometryInfo<dim>::max_children_per_face;
376 for (
unsigned int k = 0; k < n_face_points; ++k)
377 for (
unsigned int i_child = 0; i_child < this->
dofs_per_cell;
387 this->
restriction[iso][child](face * this->dofs_per_face +
390 Utilities::fixed_power<dim - 1>(.5) * q_sub.
weight(k) *
391 cached_values_face(i_child, k) *
393 face * this->dofs_per_face + i_face,
406 std::vector<AnisotropicPolynomials<dim> *> polynomials(dim);
407 for (
unsigned int dd = 0; dd < dim; ++dd)
409 std::vector<std::vector<Polynomials::Polynomial<double>>> poly(dim);
410 for (
unsigned int d = 0; d < dim; ++d)
418 const unsigned int start_cell_dofs =
425 for (
unsigned int k = 0; k < q_cell.size(); ++k)
427 for (
unsigned int d = 0; d < dim; ++d)
428 cached_values_cell(i, k, d) =
431 for (
unsigned int child = 0; child < GeometryInfo<dim>::max_children_per_cell;
436 for (
unsigned int k = 0; k < q_sub.
size(); ++k)
437 for (
unsigned int i_child = 0; i_child < this->
dofs_per_cell; ++i_child)
438 for (
unsigned int d = 0; d < dim; ++d)
439 for (
unsigned int i_weight = 0; i_weight < polynomials[d]->n();
442 this->
restriction[iso][child](start_cell_dofs + i_weight * dim +
445 q_sub.
weight(k) * cached_values_cell(i_child, k, d) *
446 polynomials[d]->compute_value(i_weight, q_sub.
point(k));
450 for (
unsigned int d = 0; d < dim; ++d)
451 delete polynomials[d];
457 std::vector<unsigned int>
463 return std::vector<unsigned int>();
472 for (
int d = 0; d < dim - 1; ++d)
473 dofs_per_face *= rt_order + 1;
476 const unsigned int interior_dofs = dim * (rt_order + 1) * dofs_per_face;
478 std::vector<unsigned int> dpo(dim + 1);
480 dpo[dim] = interior_dofs;
492 const unsigned int face_index)
const 515 return (face_index !=
537 std::vector<double> & nodal_values)
const 548 std::fill(nodal_values.begin(), nodal_values.end(), 0.);
551 for (
unsigned int face = 0; face < GeometryInfo<dim>::faces_per_cell; ++face)
552 for (
unsigned int k = 0; k < n_face_points; ++k)
557 support_point_values[face * n_face_points + k][
GeometryInfo<
558 dim>::unit_normal_direction[face]];
561 const unsigned int start_cell_dofs =
563 const unsigned int start_cell_points =
568 for (
unsigned int d = 0; d < dim; ++d)
569 nodal_values[start_cell_dofs + i * dim + d] +=
571 support_point_values[k + start_cell_points][d];
573 const unsigned int start_abf_dofs =
579 for (
unsigned int d = 0; d < dim; ++d)
580 nodal_values[start_abf_dofs + i] +=
582 support_point_values[k + start_cell_points][d];
585 for (
unsigned int face = 0; face < GeometryInfo<dim>::faces_per_cell; ++face)
589 for (
unsigned int fp = 0; fp < n_face_points; ++fp)
594 face,
false,
false,
false, n_face_points);
596 nodal_values[start_abf_dofs + i] +=
598 support_point_values[face * n_face_points + fp][
GeometryInfo<
599 dim>::unit_normal_direction[face]];
605 if (std::fabs(nodal_values[start_abf_dofs + i]) < 1.0e-16)
606 nodal_values[start_abf_dofs + i] = 0.0;
622 #include "fe_abf.inst" 624 DEAL_II_NAMESPACE_CLOSE
void initialize_restriction()
static Quadrature< dim > project_to_child(const Quadrature< dim > &quadrature, const unsigned int child_no)
std::vector< std::vector< FullMatrix< double > > > restriction
std::vector< Point< dim > > generalized_support_points
Table< 3, double > interior_weights
FullMatrix< double > interface_constraints
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
Table< 2, double > boundary_weights
virtual std::string get_name() const override
static std::vector< Polynomial< double > > generate_complete_basis(const unsigned int degree)
std::vector< Point< dim-1 > > generalized_face_support_points
static::ExceptionBase & ExcIndexRange(int arg1, int arg2, int arg3)
const Point< dim > & point(const unsigned int i) const
Table< 2, double > boundary_weights_abf
virtual std::unique_ptr< FiniteElement< dim, dim > > clone() const override
unsigned int size() const
static::ExceptionBase & ExcImpossibleInDim(int arg1)
std::vector< std::vector< FullMatrix< double > > > prolongation
void invert(const FullMatrix< number2 > &M)
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)
void reinit(const TableIndices< N > &new_size, const bool omit_default_initialization=false)
size_type size(const unsigned int i) const
Table< 3, double > interior_weights_abf
#define Assert(cond, exc)
static::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
double weight(const unsigned int i) const
virtual double shape_value_component(const unsigned int i, const Point< dim > &p, const unsigned int component) const override
static std::vector< Polynomial< number > > generate_complete_basis(const unsigned int degree)
unsigned int n_components() const
const unsigned int dofs_per_cell
static void project_to_subface(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 Quadrature< dim > project_to_all_faces(const SubQuadrature &quadrature)
static void project_to_face(const SubQuadrature &quadrature, const unsigned int face_no, std::vector< Point< dim >> &q_points)
const unsigned int rt_order
virtual std::size_t memory_consumption() const override
const unsigned int dofs_per_face
static std::vector< unsigned int > get_dpo_vector(const unsigned int degree)
void reinit_restriction_and_prolongation_matrices(const bool isotropic_restriction_only=false, const bool isotropic_prolongation_only=false)
static::ExceptionBase & ExcNotImplemented()
void initialize_support_points(const unsigned int rt_degree)
FullMatrix< double > inverse_node_matrix
virtual bool has_support_on_face(const unsigned int shape_index, const unsigned int face_index) const override
static::ExceptionBase & ExcInternalError()
static DataSetDescriptor face(const unsigned int face_no, const bool face_orientation, const bool face_flip, const bool face_rotation, const unsigned int n_quadrature_points)