17 #include <deal.II/base/std_cxx14/memory.h> 19 #include <deal.II/fe/fe_dgp_monomial.h> 20 #include <deal.II/fe/fe_tools.h> 25 DEAL_II_NAMESPACE_OPEN
48 const unsigned int start_index2d[6] = {0, 1, 4, 10, 20, 35};
49 const double points2d[35][2] = {
50 {0, 0}, {0, 0}, {1, 0}, {0, 1}, {0, 0},
51 {1, 0}, {0, 1}, {1, 1}, {0.5, 0}, {0, 0.5},
52 {0, 0}, {1, 0}, {0, 1}, {1, 1}, {1. / 3., 0},
53 {2. / 3., 0}, {0, 1. / 3.}, {0, 2. / 3.}, {0.5, 1}, {1, 0.5},
54 {0, 0}, {1, 0}, {0, 1}, {1, 1}, {0.25, 0},
55 {0.5, 0}, {0.75, 0}, {0, 0.25}, {0, 0.5}, {0, 0.75},
56 {1. / 3., 1}, {2. / 3., 1}, {1, 1. / 3.}, {1, 2. / 3.}, {0.5, 0.5}};
64 const unsigned int start_index3d[6] = {0, 1, 5, 15 };
65 const double points3d[35][3] = {{0, 0, 0},
84 generate_unit_points(
const unsigned int, std::vector<
Point<dim>> &);
88 generate_unit_points(
const unsigned int k, std::vector<
Point<1>> &p)
91 const double h = 1. / k;
92 for (
unsigned int i = 0; i < p.size(); ++i)
98 generate_unit_points(
const unsigned int k, std::vector<
Point<2>> &p)
101 Assert(p.size() == start_index2d[k + 1] - start_index2d[k],
103 for (
unsigned int i = 0; i < p.size(); ++i)
105 p[i](0) = points2d[start_index2d[k] + i][0];
106 p[i](1) = points2d[start_index2d[k] + i][1];
112 generate_unit_points(
const unsigned int k, std::vector<
Point<3>> &p)
115 Assert(p.size() == start_index3d[k + 1] - start_index3d[k],
117 for (
unsigned int i = 0; i < p.size(); ++i)
119 p[i](0) = points3d[start_index3d[k] + i][0];
120 p[i](1) = points3d[start_index3d[k] + i][1];
121 p[i](2) = points3d[start_index3d[k] + i][2];
143 std::vector<bool>(1, true)))
174 std::ostringstream namebuf;
175 namebuf <<
"FE_DGPMonomial<" << dim <<
">(" << this->
degree <<
")";
177 return namebuf.str();
183 std::unique_ptr<FiniteElement<dim, dim>>
186 return std_cxx14::make_unique<FE_DGPMonomial<dim>>(*this);
201 if (source_dgp_monomial)
206 const unsigned int m = interpolation_matrix.
m();
207 const unsigned int n = interpolation_matrix.
n();
215 const unsigned int min_mn =
216 interpolation_matrix.
m() < interpolation_matrix.
n() ?
217 interpolation_matrix.
m() :
218 interpolation_matrix.
n();
220 for (
unsigned int i = 0; i < min_mn; ++i)
221 interpolation_matrix(i, i) = 1.;
226 internal::FE_DGPMonomial::generate_unit_points(this->
degree, unit_points);
231 for (
unsigned int k = 0; k < unit_points.size(); ++k)
232 source_fe_matrix(k, j) = source_fe.
shape_value(j, unit_points[k]);
236 for (
unsigned int k = 0; k < unit_points.size(); ++k)
237 this_matrix(k, j) = this->
poly_space.compute_value(j, unit_points[k]);
241 this_matrix.
mmult(interpolation_matrix, source_fe_matrix);
261 std::vector<unsigned int>
264 std::vector<unsigned int> dpo(dim + 1, 0U);
266 for (
unsigned int i = 1; i < dim; ++i)
268 dpo[dim] *= deg + 1 + i;
287 (void)interpolation_matrix;
293 Assert(interpolation_matrix.
m() == 0,
295 Assert(interpolation_matrix.
n() == 0,
314 (void)interpolation_matrix;
320 Assert(interpolation_matrix.
m() == 0,
322 Assert(interpolation_matrix.
n() == 0,
338 std::vector<std::pair<unsigned int, unsigned int>>
345 return std::vector<std::pair<unsigned int, unsigned int>>();
349 return std::vector<std::pair<unsigned int, unsigned int>>();
356 std::vector<std::pair<unsigned int, unsigned int>>
363 return std::vector<std::pair<unsigned int, unsigned int>>();
367 return std::vector<std::pair<unsigned int, unsigned int>>();
374 std::vector<std::pair<unsigned int, unsigned int>>
381 return std::vector<std::pair<unsigned int, unsigned int>>();
385 return std::vector<std::pair<unsigned int, unsigned int>>();
412 const unsigned int face_index)
const 414 return face_index == 1 || (face_index == 0 && this->
degree == 0);
422 const unsigned int face_index)
const 424 bool support_on_face =
false;
425 if (face_index == 1 || face_index == 2)
426 support_on_face =
true;
429 const std::array<unsigned int, 2> degrees =
430 this->
poly_space.directional_degrees(shape_index);
432 if ((face_index == 0 && degrees[1] == 0) ||
433 (face_index == 3 && degrees[0] == 0))
434 support_on_face =
true;
436 return support_on_face;
444 const unsigned int face_index)
const 446 bool support_on_face =
false;
447 if (face_index == 1 || face_index == 3 || face_index == 4)
448 support_on_face =
true;
451 const std::array<unsigned int, 3> degrees =
452 this->
poly_space.directional_degrees(shape_index);
454 if ((face_index == 0 && degrees[1] == 0) ||
455 (face_index == 2 && degrees[2] == 0) ||
456 (face_index == 5 && degrees[0] == 0))
457 support_on_face =
true;
459 return support_on_face;
475 #include "fe_dgp_monomial.inst" 478 DEAL_II_NAMESPACE_CLOSE
virtual double shape_value(const unsigned int i, const Point< dim > &p) const
virtual std::vector< std::pair< unsigned int, unsigned int > > hp_quad_dof_identities(const FiniteElement< dim > &fe_other) const override
static std::vector< unsigned int > get_dpo_vector(const unsigned int degree)
virtual void get_subface_interpolation_matrix(const FiniteElement< dim > &source, const unsigned int subface, FullMatrix< double > &matrix) const override
std::vector< std::vector< FullMatrix< double > > > restriction
virtual FiniteElementDomination::Domination compare_for_face_domination(const FiniteElement< dim > &fe_other) const override
FE_DGPMonomial(const unsigned int p)
const unsigned int degree
virtual std::size_t memory_consumption() const override
virtual void get_face_interpolation_matrix(const FiniteElement< dim > &source, FullMatrix< double > &matrix) const override
#define AssertThrow(cond, exc)
virtual std::vector< std::pair< unsigned int, unsigned int > > hp_line_dof_identities(const FiniteElement< dim > &fe_other) const override
virtual void get_interpolation_matrix(const FiniteElement< dim > &source, FullMatrix< double > &matrix) const override
virtual bool has_support_on_face(const unsigned int shape_index, const unsigned int face_index) const override
void initialize_restriction()
std::vector< std::vector< FullMatrix< double > > > prolongation
#define Assert(cond, exc)
static::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
virtual std::string get_name() const =0
void mmult(FullMatrix< number2 > &C, const FullMatrix< number2 > &B, const bool adding=false) const
const unsigned int dofs_per_cell
virtual bool hp_constraints_are_implemented() const override
virtual std::vector< std::pair< unsigned int, unsigned int > > hp_vertex_dof_identities(const FiniteElement< dim > &fe_other) const override
void reinit_restriction_and_prolongation_matrices(const bool isotropic_restriction_only=false, const bool isotropic_prolongation_only=false)
static::ExceptionBase & ExcNotImplemented()
PolynomialsP< dim > poly_space
virtual std::string get_name() const override
virtual std::unique_ptr< FiniteElement< dim, dim > > clone() const override
static::ExceptionBase & ExcInternalError()