const RefinementCase<dim> &refinement_case =
RefinementCase<dim>::isotropic_refinement) const override;
+ /**
+ * @copydoc dealii::FiniteElement::get_face_interpolation_matrix()
+ */
+ void
+ get_face_interpolation_matrix(const FiniteElement<dim, spacedim> &source_fe,
+ FullMatrix<double> &interpolation_matrix,
+ const unsigned int face_no) const override;
+
+ /**
+ * @copydoc dealii::FiniteElement::get_subface_interpolation_matrix()
+ */
+ void
+ get_subface_interpolation_matrix(
+ const FiniteElement<dim, spacedim> &x_source_fe,
+ const unsigned int subface,
+ FullMatrix<double> & interpolation_matrix,
+ const unsigned int face_no) const override;
+
+ /**
+ * @copydoc dealii::FiniteElement::hp_constraints_are_implemented()
+ */
+ bool
+ hp_constraints_are_implemented() const override;
+
/**
* @copydoc dealii::FiniteElement::convert_generalized_support_point_values_to_dof_values()
*/
const unsigned int face_no,
std::vector<Point<2>> & q_points)
{
- Assert(reference_cell_type == ReferenceCell::Type::Quad, ExcNotImplemented());
- (void)reference_cell_type;
-
const unsigned int dim = 2;
AssertIndexRange(face_no, GeometryInfo<dim>::faces_per_cell);
Assert(q_points.size() == quadrature.size(),
ExcDimensionMismatch(q_points.size(), quadrature.size()));
- for (unsigned int p = 0; p < quadrature.size(); ++p)
- switch (face_no)
- {
- case 0:
- q_points[p] = Point<dim>(0, quadrature.point(p)(0));
- break;
- case 1:
- q_points[p] = Point<dim>(1, quadrature.point(p)(0));
- break;
- case 2:
- q_points[p] = Point<dim>(quadrature.point(p)(0), 0);
- break;
- case 3:
- q_points[p] = Point<dim>(quadrature.point(p)(0), 1);
- break;
- default:
- Assert(false, ExcInternalError());
- }
+ if (reference_cell_type == ReferenceCell::Type::Tri)
+ {
+ // use linear polynomial to map the reference quadrature points correctly
+ // on faces, i.e., Simplex::ScalarPolynomial<1>(1)
+ for (unsigned int p = 0; p < quadrature.size(); ++p)
+ switch (face_no)
+ {
+ case 0:
+ q_points[p] = Point<dim>(quadrature.point(p)(0), 0);
+ break;
+ case 1:
+ q_points[p] =
+ Point<dim>(1 - quadrature.point(p)(0), quadrature.point(p)(0));
+ break;
+ case 2:
+ q_points[p] = Point<dim>(0, 1 - quadrature.point(p)(0));
+ break;
+ default:
+ Assert(false, ExcInternalError());
+ }
+ }
+ else if (reference_cell_type == ReferenceCell::Type::Quad)
+ {
+ for (unsigned int p = 0; p < quadrature.size(); ++p)
+ switch (face_no)
+ {
+ case 0:
+ q_points[p] = Point<dim>(0, quadrature.point(p)(0));
+ break;
+ case 1:
+ q_points[p] = Point<dim>(1, quadrature.point(p)(0));
+ break;
+ case 2:
+ q_points[p] = Point<dim>(quadrature.point(p)(0), 0);
+ break;
+ case 3:
+ q_points[p] = Point<dim>(quadrature.point(p)(0), 1);
+ break;
+ default:
+ Assert(false, ExcInternalError());
+ }
+ }
+ else
+ {
+ Assert(false, ExcInternalError());
+ }
}
std::vector<Point<2>> & q_points,
const RefinementCase<1> &)
{
- Assert(reference_cell_type == ReferenceCell::Type::Quad, ExcNotImplemented());
- (void)reference_cell_type;
-
const unsigned int dim = 2;
AssertIndexRange(face_no, GeometryInfo<dim>::faces_per_cell);
AssertIndexRange(subface_no, GeometryInfo<dim>::max_children_per_face);
Assert(q_points.size() == quadrature.size(),
ExcDimensionMismatch(q_points.size(), quadrature.size()));
- for (unsigned int p = 0; p < quadrature.size(); ++p)
- switch (face_no)
- {
- case 0:
- switch (subface_no)
- {
- case 0:
- q_points[p] = Point<dim>(0, quadrature.point(p)(0) / 2);
- break;
- case 1:
- q_points[p] = Point<dim>(0, quadrature.point(p)(0) / 2 + 0.5);
- break;
- default:
- Assert(false, ExcInternalError());
- }
- break;
- case 1:
- switch (subface_no)
- {
- case 0:
- q_points[p] = Point<dim>(1, quadrature.point(p)(0) / 2);
- break;
- case 1:
- q_points[p] = Point<dim>(1, quadrature.point(p)(0) / 2 + 0.5);
- break;
- default:
- Assert(false, ExcInternalError());
- }
- break;
- case 2:
- switch (subface_no)
- {
- case 0:
- q_points[p] = Point<dim>(quadrature.point(p)(0) / 2, 0);
- break;
- case 1:
- q_points[p] = Point<dim>(quadrature.point(p)(0) / 2 + 0.5, 0);
- break;
- default:
- Assert(false, ExcInternalError());
- }
- break;
- case 3:
- switch (subface_no)
- {
- case 0:
- q_points[p] = Point<dim>(quadrature.point(p)(0) / 2, 1);
- break;
- case 1:
- q_points[p] = Point<dim>(quadrature.point(p)(0) / 2 + 0.5, 1);
- break;
- default:
- Assert(false, ExcInternalError());
- }
- break;
-
- default:
- Assert(false, ExcInternalError());
- }
+ if (reference_cell_type == ReferenceCell::Type::Tri)
+ {
+ // use linear polynomial to map the reference quadrature points correctly
+ // on faces, i.e., Simplex::ScalarPolynomial<1>(1)
+ for (unsigned int p = 0; p < quadrature.size(); ++p)
+ switch (face_no)
+ {
+ case 0:
+ switch (subface_no)
+ {
+ case 0:
+ q_points[p] = Point<dim>(quadrature.point(p)(0) / 2, 0);
+ break;
+ case 1:
+ q_points[p] =
+ Point<dim>(0.5 + quadrature.point(p)(0) / 2, 0);
+ break;
+ default:
+ Assert(false, ExcInternalError());
+ }
+ break;
+ case 1:
+ switch (subface_no)
+ {
+ case 0:
+ q_points[p] = Point<dim>(1 - quadrature.point(p)(0) / 2,
+ quadrature.point(p)(0) / 2);
+ break;
+ case 1:
+ q_points[p] = Point<dim>(0.5 - quadrature.point(p)(0) / 2,
+ 0.5 + quadrature.point(p)(0) / 2);
+ break;
+ default:
+ Assert(false, ExcInternalError());
+ }
+ break;
+ case 2:
+ switch (subface_no)
+ {
+ case 0:
+ q_points[p] = Point<dim>(0, 1 - quadrature.point(p)(0) / 2);
+ break;
+ case 1:
+ q_points[p] =
+ Point<dim>(0, 0.5 - quadrature.point(p)(0) / 2);
+ break;
+ default:
+ Assert(false, ExcInternalError());
+ }
+ break;
+ default:
+ Assert(false, ExcInternalError());
+ }
+ }
+ else if (reference_cell_type == ReferenceCell::Type::Quad)
+ {
+ for (unsigned int p = 0; p < quadrature.size(); ++p)
+ switch (face_no)
+ {
+ case 0:
+ switch (subface_no)
+ {
+ case 0:
+ q_points[p] = Point<dim>(0, quadrature.point(p)(0) / 2);
+ break;
+ case 1:
+ q_points[p] =
+ Point<dim>(0, quadrature.point(p)(0) / 2 + 0.5);
+ break;
+ default:
+ Assert(false, ExcInternalError());
+ }
+ break;
+ case 1:
+ switch (subface_no)
+ {
+ case 0:
+ q_points[p] = Point<dim>(1, quadrature.point(p)(0) / 2);
+ break;
+ case 1:
+ q_points[p] =
+ Point<dim>(1, quadrature.point(p)(0) / 2 + 0.5);
+ break;
+ default:
+ Assert(false, ExcInternalError());
+ }
+ break;
+ case 2:
+ switch (subface_no)
+ {
+ case 0:
+ q_points[p] = Point<dim>(quadrature.point(p)(0) / 2, 0);
+ break;
+ case 1:
+ q_points[p] =
+ Point<dim>(quadrature.point(p)(0) / 2 + 0.5, 0);
+ break;
+ default:
+ Assert(false, ExcInternalError());
+ }
+ break;
+ case 3:
+ switch (subface_no)
+ {
+ case 0:
+ q_points[p] = Point<dim>(quadrature.point(p)(0) / 2, 1);
+ break;
+ case 1:
+ q_points[p] =
+ Point<dim>(quadrature.point(p)(0) / 2 + 0.5, 1);
+ break;
+ default:
+ Assert(false, ExcInternalError());
+ }
+ break;
+
+ default:
+ Assert(false, ExcInternalError());
+ }
+ }
+ else
+ {
+ Assert(false, ExcInternalError());
+ }
}
{
if (reference_cell_type == ReferenceCell::Type::Tri)
{
+ const auto support_points_line =
+ [](const auto &face, const auto &orientation) -> std::vector<Point<2>> {
+ std::array<Point<2>, 2> vertices;
+ std::copy_n(face.first.begin(), face.first.size(), vertices.begin());
+ const auto temp =
+ ReferenceCell::Type(ReferenceCell::Type::Line)
+ .permute_according_orientation(vertices, orientation);
+ return std::vector<Point<2>>(temp.begin(),
+ temp.begin() + face.first.size());
+ };
+
// reference faces (defined by its support points and arc length)
const std::array<std::pair<std::array<Point<2>, 2>, double>, 3> faces = {
{{{{Point<2>(0.0, 0.0), Point<2>(1.0, 0.0)}}, 1.0},
{
const auto &face = faces[face_no];
- std::array<Point<2>, 2> support_points;
-
// determine support point of the current line with the correct
// orientation
- switch (orientation)
- {
- case 0:
- support_points = {{face.first[1], face.first[0]}};
- break;
- case 1:
- support_points = {{face.first[0], face.first[1]}};
- break;
- default:
- Assert(false, ExcNotImplemented());
- }
+ std::vector<Point<2>> support_points =
+ support_points_line(face, orientation);
// the quadrature rule to be projected ...
const auto &sub_quadrature_points =
}
}
+#ifdef DEBUG
// make sure that the row sum of each of the matrices is 1 at this
// point. this must be so since the shape functions sum up to 1
for (unsigned int j = 0; j < source_fe->n_dofs_per_face(face_no); ++j)
Assert(std::fabs(sum - 1) < eps, ExcInternalError());
}
+#endif
}
else if (dynamic_cast<const FE_Nothing<dim> *>(&x_source_fe) != nullptr)
{
#include <deal.II/base/config.h>
+#include <deal.II/base/qprojector.h>
+
#include <deal.II/fe/fe_dgq.h>
#include <deal.II/fe/fe_nothing.h>
#include <deal.II/fe/fe_q.h>
+ template <int dim, int spacedim>
+ void
+ FE_Poly<dim, spacedim>::get_face_interpolation_matrix(
+ const FiniteElement<dim, spacedim> &x_source_fe,
+ FullMatrix<double> & interpolation_matrix,
+ const unsigned int face_no) const
+ {
+ Assert(interpolation_matrix.m() == x_source_fe.n_dofs_per_face(face_no),
+ ExcDimensionMismatch(interpolation_matrix.m(),
+ x_source_fe.n_dofs_per_face(face_no)));
+
+ if (const FE_Poly<dim, spacedim> *source_fe =
+ dynamic_cast<const FE_Poly<dim, spacedim> *>(&x_source_fe))
+ {
+ const Quadrature<dim - 1> quad_face_support(
+ source_fe->get_unit_face_support_points(face_no));
+
+ const double eps = 2e-13 * this->degree * (dim - 1);
+
+ std::vector<Point<dim>> face_quadrature_points(
+ quad_face_support.size());
+ QProjector<dim>::project_to_face(this->reference_cell_type(),
+ quad_face_support,
+ face_no,
+ face_quadrature_points);
+
+ for (unsigned int i = 0; i < source_fe->n_dofs_per_face(face_no); ++i)
+ for (unsigned int j = 0; j < this->n_dofs_per_face(face_no); ++j)
+ {
+ double matrix_entry =
+ this->shape_value(this->face_to_cell_index(j, 0),
+ face_quadrature_points[i]);
+
+ // Correct the interpolated value. I.e. if it is close to 1 or
+ // 0, make it exactly 1 or 0. Unfortunately, this is required to
+ // avoid problems with higher order elements.
+ if (std::fabs(matrix_entry - 1.0) < eps)
+ matrix_entry = 1.0;
+ if (std::fabs(matrix_entry) < eps)
+ matrix_entry = 0.0;
+
+ interpolation_matrix(i, j) = matrix_entry;
+ }
+
+#ifdef DEBUG
+ for (unsigned int j = 0; j < source_fe->n_dofs_per_face(face_no); ++j)
+ {
+ double sum = 0.;
+
+ for (unsigned int i = 0; i < this->n_dofs_per_face(face_no); ++i)
+ sum += interpolation_matrix(j, i);
+
+ Assert(std::fabs(sum - 1) < eps, ExcInternalError());
+ }
+#endif
+ }
+ else if (dynamic_cast<const FE_Nothing<dim> *>(&x_source_fe) != nullptr)
+ {
+ // nothing to do here, the FE_Nothing has no degrees of freedom anyway
+ }
+ else
+ AssertThrow(
+ false,
+ (typename FiniteElement<dim,
+ spacedim>::ExcInterpolationNotImplemented()));
+ }
+
+
+
+ template <int dim, int spacedim>
+ void
+ FE_Poly<dim, spacedim>::get_subface_interpolation_matrix(
+ const FiniteElement<dim, spacedim> &x_source_fe,
+ const unsigned int subface,
+ FullMatrix<double> & interpolation_matrix,
+ const unsigned int face_no) const
+ {
+ Assert(interpolation_matrix.m() == x_source_fe.n_dofs_per_face(face_no),
+ ExcDimensionMismatch(interpolation_matrix.m(),
+ x_source_fe.n_dofs_per_face(face_no)));
+
+ if (const FE_Poly<dim, spacedim> *source_fe =
+ dynamic_cast<const FE_Poly<dim, spacedim> *>(&x_source_fe))
+ {
+ const Quadrature<dim - 1> quad_face_support(
+ source_fe->get_unit_face_support_points(face_no));
+
+ const double eps = 2e-13 * this->degree * (dim - 1);
+
+ std::vector<Point<dim>> subface_quadrature_points(
+ quad_face_support.size());
+ QProjector<dim>::project_to_subface(this->reference_cell_type(),
+ quad_face_support,
+ face_no,
+ subface,
+ subface_quadrature_points);
+
+ for (unsigned int i = 0; i < source_fe->n_dofs_per_face(face_no); ++i)
+ for (unsigned int j = 0; j < this->n_dofs_per_face(face_no); ++j)
+ {
+ double matrix_entry =
+ this->shape_value(this->face_to_cell_index(j, 0),
+ subface_quadrature_points[i]);
+
+ // Correct the interpolated value. I.e. if it is close to 1 or
+ // 0, make it exactly 1 or 0. Unfortunately, this is required to
+ // avoid problems with higher order elements.
+ if (std::fabs(matrix_entry - 1.0) < eps)
+ matrix_entry = 1.0;
+ if (std::fabs(matrix_entry) < eps)
+ matrix_entry = 0.0;
+
+ interpolation_matrix(i, j) = matrix_entry;
+ }
+
+#ifdef DEBUG
+ for (unsigned int j = 0; j < source_fe->n_dofs_per_face(face_no); ++j)
+ {
+ double sum = 0.;
+
+ for (unsigned int i = 0; i < this->n_dofs_per_face(face_no); ++i)
+ sum += interpolation_matrix(j, i);
+
+ Assert(std::fabs(sum - 1) < eps, ExcInternalError());
+ }
+#endif
+ }
+ else if (dynamic_cast<const FE_Nothing<dim> *>(&x_source_fe) != nullptr)
+ {
+ // nothing to do here, the FE_Nothing has no degrees of freedom anyway
+ }
+ else
+ AssertThrow(
+ false,
+ (typename FiniteElement<dim,
+ spacedim>::ExcInterpolationNotImplemented()));
+ }
+
+
+
+ template <int dim, int spacedim>
+ bool
+ FE_Poly<dim, spacedim>::hp_constraints_are_implemented() const
+ {
+ return true;
+ }
+
+
+
template <int dim, int spacedim>
void
FE_Poly<dim, spacedim>::
// verify hanging node constraints on locally h-refined simplex mesh
+//
+// dofs will be enumerated as follows
+// scenario 1:
+// 1-------0
+// |\ |
+// | \ |
+// 5---6 |
+// |\ |\ |
+// | \| \|
+// 3---4---2
+
#include <deal.II/dofs/dof_handler.h>
#include <deal.II/dofs/dof_tools.h>
+#include <deal.II/fe/mapping_fe.h>
+
#include <deal.II/grid/grid_out.h>
#include <deal.II/grid/tria.h>
#include "../tests.h"
+// ----- diagnostics -----
+
+template <int dim>
+void
+print_dof_indices_on_faces(const DoFHandler<dim> &dofh)
+{
+ std::vector<types::global_dof_index> dof_indices;
+
+ for (const auto &cell : dofh.active_cell_iterators())
+ for (unsigned int f = 0; f < cell->n_faces(); ++f)
+ {
+ const auto &face = cell->face(f);
+
+ if (face->has_children())
+ {
+ for (unsigned int sf = 0; sf < face->n_children(); ++sf)
+ {
+ const auto &subface = face->child(sf);
+ Assert(subface->n_active_fe_indices() == 1, ExcInternalError());
+ const unsigned int subface_fe_index =
+ subface->nth_active_fe_index(0);
+ const auto &subface_fe = dofh.get_fe(subface_fe_index);
+
+ dof_indices.resize(subface_fe.n_dofs_per_face(f));
+ subface->get_dof_indices(dof_indices, subface_fe_index);
+
+ deallog << "cell:" << cell->active_cell_index() << " face:" << f
+ << " subface:" << sf << " dofs:";
+ for (const auto &i : dof_indices)
+ deallog << i << " ";
+ deallog << std::endl;
+ }
+ }
+ else
+ {
+ Assert(face->n_active_fe_indices() == 1, ExcInternalError());
+ const unsigned int face_fe_index = face->nth_active_fe_index(0);
+ const auto & face_fe = dofh.get_fe(face_fe_index);
+
+ dof_indices.resize(face_fe.n_dofs_per_face(f));
+ face->get_dof_indices(dof_indices, face_fe_index);
+
+ deallog << "cell:" << cell->active_cell_index() << " face:" << f
+ << " dofs:";
+ for (const auto &i : dof_indices)
+ deallog << i << " ";
+ deallog << std::endl;
+ }
+ }
+}
+
+
+template <int dim>
+void
+print_dof_points(const DoFHandler<dim> &dofh)
+{
+ std::vector<Point<dim>> points(dofh.n_dofs());
+ DoFTools::map_dofs_to_support_points(MappingFE<dim>(dofh.get_fe()),
+ dofh,
+ points);
+
+ for (unsigned int i = 0; i < dofh.n_dofs(); ++i)
+ deallog << "dof:" << i << " point:" << points[i] << std::endl;
+}
+
+
+
+// ----- test -----
+
template <int dim>
void
test()
dofh.distribute_dofs(Simplex::FE_P<dim>(1));
deallog << "ndofs: " << dofh.n_dofs() << std::endl;
+#if false
+ print_dof_points(dofh);
+ print_dof_indices_on_faces(dofh);
+#endif
+
// hanging node constraints
AffineConstraints<double> constraints;
- // DoFTools::make_hanging_node_constraints(dofh, constraints);
+ DoFTools::make_hanging_node_constraints(dofh, constraints);
constraints.print(deallog.get_file_stream());
deallog << "OK" << std::endl;
DEAL:2d::ndofs: 7
+ 6 2: 0.500000
+ 6 1: 0.500000
DEAL:2d::OK
// verify hanging node constraints on locally p-refined simplex mesh
+//
+// dofs will be enumerated as follows
+// scenario 1: scenario 2:
+// 6-------4 2---4---3
+// |\ | |\ |
+// | \ | | \ |
+// 3 2 | | 5 6
+// | \ | | \ |
+// | \| | \|
+// 0---1---5 0-------1
+
#include <deal.II/dofs/dof_handler.h>
#include <deal.II/dofs/dof_tools.h>
+#include <deal.II/fe/mapping_fe.h>
+
#include <deal.II/grid/tria.h>
#include <deal.II/hp/fe_collection.h>
+#include <deal.II/hp/mapping_collection.h>
#include <deal.II/lac/affine_constraints.h>
#include "../tests.h"
+// ----- diagnostics -----
+
+template <int dim>
+void
+print_dof_indices_on_faces(const DoFHandler<dim> &dofh)
+{
+ std::vector<types::global_dof_index> dof_indices;
+
+ for (const auto &cell : dofh.active_cell_iterators())
+ for (unsigned int f = 0; f < cell->n_faces(); ++f)
+ {
+ const auto &face = cell->face(f);
+
+ Assert(!face->has_children(), ExcInternalError());
+
+ const unsigned int fe_index = cell->active_fe_index();
+ const auto & fe = cell->get_fe();
+
+ dof_indices.resize(fe.n_dofs_per_face(f));
+ face->get_dof_indices(dof_indices, fe_index);
+
+ deallog << "cell:" << cell->active_cell_index() << " face:" << f
+ << " dofs:";
+ for (const auto &i : dof_indices)
+ deallog << i << " ";
+ deallog << std::endl;
+ }
+}
+
+
+template <int dim>
+void
+print_dof_points(const DoFHandler<dim> &dofh)
+{
+ hp::MappingCollection<dim> mapping;
+ for (unsigned int i = 0; i < dofh.get_fe_collection().size(); ++i)
+ mapping.push_back(MappingFE<dim>(dofh.get_fe(i)));
+
+ std::vector<Point<dim>> points(dofh.n_dofs());
+ DoFTools::map_dofs_to_support_points(mapping, dofh, points);
+
+ for (unsigned int i = 0; i < dofh.n_dofs(); ++i)
+ deallog << "dof:" << i << " point:" << points[i] << std::endl;
+}
+
+
+// ----- test -----
+
template <int dim>
void
test(const hp::FECollection<dim> &fes)
Triangulation<dim> tria;
GridGenerator::subdivided_hyper_cube_with_simplices(tria, 1);
+#if false
+ GridOut grid_out;
+ grid_out.write_vtk(tria, deallog.get_file_stream());
+#endif
+
DoFHandler<dim> dofh(tria);
dofh.begin_active()->set_active_fe_index(1);
dofh.distribute_dofs(fes);
deallog << "ndofs: " << dofh.n_dofs() << std::endl;
+#if false
+ print_dof_points(dofh);
+ print_dof_indices_on_faces(dofh);
+#endif
+
// hanging node constraints
AffineConstraints<double> constraints;
- // DoFTools::make_hanging_node_constraints(dofh, constraints);
+ DoFTools::make_hanging_node_constraints(dofh, constraints);
constraints.print(deallog.get_file_stream());
deallog << "OK" << std::endl;
DEAL:2d::ndofs: 7
+ 2 6: 0.500000
+ 2 5: 0.500000
DEAL:2d::OK
DEAL:2d::ndofs: 7
+ 5 1: 0.500000
+ 5 2: 0.500000
DEAL:2d::OK
projects / dealii.git / commitdiff