standard_line_to_face_and_line_index(
const unsigned int line) const override
{
- Assert(false, ExcNotImplemented());
-
static const std::array<unsigned int, 2> table[6] = {
{{0, 0}}, {{0, 1}}, {{0, 2}}, {{1, 1}}, {{1, 2}}, {{2, 1}}};
standard_line_to_face_and_line_index(
const unsigned int line) const override
{
- Assert(false, ExcNotImplemented());
-
- static const std::array<unsigned int, 2> table[6] = {
- {{0, 0}}, {{0, 1}}, {{0, 2}}, {{1, 1}}, {{1, 2}}, {{2, 1}}};
+ static const std::array<unsigned int, 2> table[9] = {{{0, 0}},
+ {{0, 2}},
+ {{0, 1}},
+ {{1, 0}},
+ {{1, 1}},
+ {{1, 2}},
+ {{2, 0}},
+ {{2, 1}},
+ {{3, 1}}};
return table[line];
}
const unsigned int face,
const unsigned char face_orientation) const override
{
- Assert(false, ExcNotImplemented());
-
- (void)face;
-
- static const std::array<std::array<unsigned int, 3>, 6> table = {
- {{{2, 1, 0}},
- {{0, 1, 2}},
- {{1, 2, 0}},
- {{0, 2, 1}},
- {{1, 0, 2}},
- {{2, 0, 1}}}};
+ if (face > 1) // QUAD
+ {
+ return GeometryInfo<3>::standard_to_real_face_line(
+ line,
+ get_bit(face_orientation, 0),
+ get_bit(face_orientation, 2),
+ get_bit(face_orientation, 1));
+ }
+ else // TRI
+ {
+ static const std::array<std::array<unsigned int, 3>, 6> table = {
+ {{{2, 1, 0}},
+ {{0, 1, 2}},
+ {{1, 2, 0}},
+ {{0, 2, 1}},
+ {{1, 0, 2}},
+ {{2, 0, 1}}}};
- return table[face_orientation][line];
+ return table[face_orientation][line];
+ }
}
bool
const FiniteElement<dim, spacedim> &fe_other) const override;
};
+
+ /**
+ * Implementation of a scalar Lagrange finite element on a wedge that yields
+ * the finite element space of continuous, piecewise polynomials of
+ * degree p.
+ *
+ * @ingroup simplex
+ */
+ template <int dim, int spacedim = dim>
+ class FE_WedgeP : public dealii::FE_Poly<dim, spacedim>
+ {
+ public:
+ /**
+ * Constructor.
+ */
+ FE_WedgeP(const unsigned int degree);
+
+ /**
+ * @copydoc dealii::FiniteElement::clone()
+ */
+ std::unique_ptr<FiniteElement<dim, spacedim>>
+ clone() const override;
+
+ /**
+ * Return a string that uniquely identifies a finite element. This class
+ * returns <tt>Simplex::FE_WedgeP<dim>(degree)</tt>, with @p dim and @p degree
+ * replaced by appropriate values.
+ */
+ std::string
+ get_name() const override;
+
+ /**
+ * @copydoc dealii::FiniteElement::compare_for_domination()
+ */
+ FiniteElementDomination::Domination
+ compare_for_domination(const FiniteElement<dim, spacedim> &fe_other,
+ const unsigned int codim) const override;
+
+ /**
+ * @copydoc dealii::FiniteElement::hp_vertex_dof_identities()
+ */
+ std::vector<std::pair<unsigned int, unsigned int>>
+ hp_vertex_dof_identities(
+ const FiniteElement<dim, spacedim> &fe_other) const override;
+
+ /**
+ * @copydoc dealii::FiniteElement::hp_line_dof_identities()
+ */
+ std::vector<std::pair<unsigned int, unsigned int>>
+ hp_line_dof_identities(
+ const FiniteElement<dim, spacedim> &fe_other) const override;
+
+ /**
+ * @copydoc dealii::FiniteElement::hp_quad_dof_identities()
+ */
+ std::vector<std::pair<unsigned int, unsigned int>>
+ hp_quad_dof_identities(const FiniteElement<dim, spacedim> &fe_other,
+ const unsigned int face_no = 0) const override;
+ };
+
+
} // namespace Simplex
DEAL_II_NAMESPACE_CLOSE
clone() const override;
};
+
+
+ /**
+ * Polynomials defined on wedge entities. This class is basis of
+ * Simplex::FE_WedgeP.
+ *
+ * The polynomials are created via a tensor product of a
+ * Simplex::ScalarPolynomial<2>(degree) and a
+ * Simplex::ScalarPolynomial<1>(degree), however, are re-numerated to better
+ * match the definition of FiniteElement.
+ */
+ template <int dim>
+ class ScalarWedgePolynomial : public ScalarPolynomialsBase<dim>
+ {
+ public:
+ /**
+ * Make the dimension available to the outside.
+ */
+ static const unsigned int dimension = dim;
+
+ /*
+ * Constructor taking the polynomial @p degree as input.
+ *
+ * @note Currently, only linear (degree=1) and quadratic polynomials
+ * (degree=2) are implemented.
+ */
+ ScalarWedgePolynomial(const unsigned int degree);
+
+ /**
+ * @copydoc ScalarPolynomialsBase::evaluate()
+ *
+ * @note Currently, only the vectors @p values and @p grads are filled.
+ */
+ void
+ evaluate(const Point<dim> & unit_point,
+ std::vector<double> & values,
+ std::vector<Tensor<1, dim>> &grads,
+ std::vector<Tensor<2, dim>> &grad_grads,
+ std::vector<Tensor<3, dim>> &third_derivatives,
+ std::vector<Tensor<4, dim>> &fourth_derivatives) const override;
+
+ /**
+ * @copydoc ScalarPolynomialsBase::compute_value()
+ */
+ double
+ compute_value(const unsigned int i, const Point<dim> &p) const override;
+
+ /**
+ * @copydoc ScalarPolynomialsBase::compute_derivative()
+ *
+ * @note Currently, only implemented for first derivative.
+ */
+ template <int order>
+ Tensor<order, dim>
+ compute_derivative(const unsigned int i, const Point<dim> &p) const;
+
+ /**
+ * @copydoc ScalarPolynomialsBase::compute_1st_derivative()
+ */
+ Tensor<1, dim>
+ compute_1st_derivative(const unsigned int i,
+ const Point<dim> & p) const override;
+
+ /**
+ * @copydoc ScalarPolynomialsBase::compute_2nd_derivative()
+ *
+ * @note Not implemented yet.
+ */
+ Tensor<2, dim>
+ compute_2nd_derivative(const unsigned int i,
+ const Point<dim> & p) const override;
+
+ /**
+ * @copydoc ScalarPolynomialsBase::compute_3rd_derivative()
+ *
+ * @note Not implemented yet.
+ */
+ Tensor<3, dim>
+ compute_3rd_derivative(const unsigned int i,
+ const Point<dim> & p) const override;
+
+ /**
+ * @copydoc ScalarPolynomialsBase::compute_4th_derivative()
+ *
+ * @note Not implemented yet.
+ */
+ Tensor<4, dim>
+ compute_4th_derivative(const unsigned int i,
+ const Point<dim> & p) const override;
+
+ /**
+ * @copydoc ScalarPolynomialsBase::compute_grad()
+ *
+ * @note Not implemented yet.
+ */
+ Tensor<1, dim>
+ compute_grad(const unsigned int i, const Point<dim> &p) const override;
+
+ /**
+ * @copydoc ScalarPolynomialsBase::compute_grad_grad()
+ *
+ * @note Not implemented yet.
+ */
+ Tensor<2, dim>
+ compute_grad_grad(const unsigned int i, const Point<dim> &p) const override;
+
+ /**
+ * @copydoc ScalarPolynomialsBase::name()
+ */
+ std::string
+ name() const override;
+
+ /**
+ * @copydoc ScalarPolynomialsBase::clone()
+ */
+ virtual std::unique_ptr<ScalarPolynomialsBase<dim>>
+ clone() const override;
+ };
+
+
+
// template functions
template <int dim>
template <int order>
return der;
}
+
+
+
+ template <int dim>
+ template <int order>
+ Tensor<order, dim>
+ ScalarWedgePolynomial<dim>::compute_derivative(const unsigned int i,
+ const Point<dim> & p) const
+ {
+ Tensor<order, dim> der;
+
+ AssertDimension(order, 1);
+ const auto grad = compute_grad(i, p);
+
+ for (unsigned int i = 0; i < dim; i++)
+ der[i] = grad[i];
+
+ return der;
+ }
+
} // namespace Simplex
DEAL_II_NAMESPACE_CLOSE
return dpo;
}
+
+ /**
+ * Helper function to set up the dpo vector of FE_WedgeP for a given @p degree.
+ */
+ internal::GenericDoFsPerObject
+ get_dpo_vector_fe_wedge(const unsigned int degree)
+ {
+ internal::GenericDoFsPerObject dpo;
+
+ // all dofs are internal
+ if (degree == 1)
+ {
+ dpo.dofs_per_object_exclusive = {{1}, {0}, {0, 0, 0, 0, 0}, {0}};
+ dpo.dofs_per_object_inclusive = {{1}, {2}, {3, 3, 4, 4, 4}, {6}};
+ dpo.object_index = {{}, {6}, {6}, {6}};
+ dpo.first_object_index_on_face = {{},
+ {3, 3, 4, 4, 4},
+ {3, 3, 4, 4, 4}};
+ }
+ else if (degree == 2)
+ {
+ dpo.dofs_per_object_exclusive = {{1}, {1}, {0, 0, 1, 1, 1}, {0}};
+ dpo.dofs_per_object_inclusive = {{1}, {3}, {6, 6, 9, 9, 9}, {18}};
+ dpo.object_index = {{}, {6}, {15, 15, 15, 16, 17}, {18}};
+ dpo.first_object_index_on_face = {{},
+ {3, 3, 4, 4, 4},
+ {6, 6, 8, 8, 8}};
+ }
+ else
+ {
+ Assert(false, ExcNotImplemented());
+ }
+
+ return dpo;
+ }
} // namespace
+ template <int dim, int spacedim>
+ FE_WedgeP<dim, spacedim>::FE_WedgeP(const unsigned int degree)
+ : dealii::FE_Poly<dim, spacedim>(
+ Simplex::ScalarWedgePolynomial<dim>(degree),
+ FiniteElementData<dim>(get_dpo_vector_fe_wedge(degree),
+ ReferenceCell::Type::Wedge,
+ 1,
+ degree,
+ FiniteElementData<dim>::L2),
+ std::vector<bool>(FiniteElementData<dim>(get_dpo_vector_fe_wedge(
+ degree),
+ ReferenceCell::Type::Wedge,
+ 1,
+ degree)
+ .dofs_per_cell,
+ true),
+ std::vector<ComponentMask>(
+ FiniteElementData<dim>(get_dpo_vector_fe_wedge(degree),
+ ReferenceCell::Type::Wedge,
+ 1,
+ degree)
+ .dofs_per_cell,
+ std::vector<bool>(1, true)))
+ {
+ AssertDimension(dim, 3);
+
+
+ if (degree == 1)
+ {
+ this->unit_support_points.emplace_back(1.0, 0.0, 0.0);
+ this->unit_support_points.emplace_back(0.0, 1.0, 0.0);
+ this->unit_support_points.emplace_back(0.0, 0.0, 0.0);
+ this->unit_support_points.emplace_back(1.0, 0.0, 1.0);
+ this->unit_support_points.emplace_back(0.0, 1.0, 1.0);
+ this->unit_support_points.emplace_back(0.0, 0.0, 1.0);
+
+ // TODO
+ this->unit_face_support_points[0].emplace_back(1.0, 0.0);
+ this->unit_face_support_points[0].emplace_back(0.0, 1.0);
+ this->unit_face_support_points[0].emplace_back(0.0, 0.0);
+ }
+ }
+
+
+
+ template <int dim, int spacedim>
+ std::unique_ptr<FiniteElement<dim, spacedim>>
+ FE_WedgeP<dim, spacedim>::clone() const
+ {
+ return std::make_unique<FE_WedgeP<dim, spacedim>>(*this);
+ }
+
+
+
+ template <int dim, int spacedim>
+ std::string
+ FE_WedgeP<dim, spacedim>::get_name() const
+ {
+ std::ostringstream namebuf;
+ namebuf << "FE_WedgeP<" << dim << ">(" << this->degree << ")";
+
+ return namebuf.str();
+ }
+
+
+
+ template <int dim, int spacedim>
+ FiniteElementDomination::Domination
+ FE_WedgeP<dim, spacedim>::compare_for_domination(
+ const FiniteElement<dim, spacedim> &fe_other,
+ const unsigned int codim) const
+ {
+ (void)fe_other;
+ (void)codim;
+
+ Assert((dynamic_cast<const Simplex::FE_P<dim, spacedim> *>(&fe_other)) ||
+ (dynamic_cast<const FE_Q<dim, spacedim> *>(&fe_other)),
+ ExcNotImplemented());
+
+ return FiniteElementDomination::this_element_dominates;
+ }
+
+
+
+ template <int dim, int spacedim>
+ std::vector<std::pair<unsigned int, unsigned int>>
+ FE_WedgeP<dim, spacedim>::hp_vertex_dof_identities(
+ const FiniteElement<dim, spacedim> &fe_other) const
+ {
+ (void)fe_other;
+
+ Assert((dynamic_cast<const Simplex::FE_P<dim, spacedim> *>(&fe_other)) ||
+ (dynamic_cast<const FE_Q<dim, spacedim> *>(&fe_other)),
+ ExcNotImplemented());
+
+ return {{0, 0}};
+ }
+
+
+
+ template <int dim, int spacedim>
+ std::vector<std::pair<unsigned int, unsigned int>>
+ FE_WedgeP<dim, spacedim>::hp_line_dof_identities(
+ const FiniteElement<dim, spacedim> &fe_other) const
+ {
+ (void)fe_other;
+
+ Assert((dynamic_cast<const Simplex::FE_P<dim, spacedim> *>(&fe_other)) ||
+ (dynamic_cast<const FE_Q<dim, spacedim> *>(&fe_other)),
+ ExcNotImplemented());
+
+ std::vector<std::pair<unsigned int, unsigned int>> result;
+
+ for (unsigned int i = 0; i < this->degree - 1; ++i)
+ result.emplace_back(i, i);
+
+ return result;
+ }
+
+
+
+ template <int dim, int spacedim>
+ std::vector<std::pair<unsigned int, unsigned int>>
+ FE_WedgeP<dim, spacedim>::hp_quad_dof_identities(
+ const FiniteElement<dim, spacedim> &fe_other,
+ const unsigned int face_no) const
+ {
+ (void)fe_other;
+
+
+ AssertIndexRange(face_no, 5);
+
+ if (face_no < 2)
+ {
+ Assert((dynamic_cast<const Simplex::FE_P<dim, spacedim> *>(&fe_other)),
+ ExcNotImplemented());
+ }
+ else
+ {
+ Assert((dynamic_cast<const FE_Q<dim, spacedim> *>(&fe_other)),
+ ExcNotImplemented());
+ }
+
+ std::vector<std::pair<unsigned int, unsigned int>> result;
+
+ for (unsigned int i = 0; i < this->n_dofs_per_quad(face_no); ++i)
+ result.emplace_back(i, i);
+
+ return result;
+ }
+
+
+
} // namespace Simplex
// explicit instantiations
template class Simplex::FE_Poly<deal_II_dimension, deal_II_space_dimension>;
template class Simplex::FE_P<deal_II_dimension, deal_II_space_dimension>;
template class Simplex::FE_DGP<deal_II_dimension, deal_II_space_dimension>;
+ template class Simplex::FE_WedgeP<deal_II_dimension,
+ deal_II_space_dimension>;
#endif
}
return 0;
}
+
+ unsigned int
+ compute_n_polynomials_wedge(const unsigned int dim,
+ const unsigned int degree)
+ {
+ if (dim == 3)
+ {
+ if (degree == 1)
+ return 6;
+ if (degree == 2)
+ return 18;
+ }
+
+ Assert(false, ExcNotImplemented());
+
+ return 0;
+ }
} // namespace
return std::make_unique<ScalarPolynomial<dim>>(*this);
}
+
+
+ template <int dim>
+ ScalarWedgePolynomial<dim>::ScalarWedgePolynomial(const unsigned int degree)
+ : ScalarPolynomialsBase<dim>(degree,
+ compute_n_polynomials_wedge(dim, degree))
+ {}
+
+
+ namespace
+ {
+ /**
+ * TODO
+ */
+ static const constexpr std::array<std::array<unsigned int, 2>, 6>
+ wedge_table_1{
+ {{{0, 0}}, {{1, 0}}, {{2, 0}}, {{0, 1}}, {{1, 1}}, {{2, 1}}}};
+
+ /**
+ * TODO
+ */
+ static const constexpr std::array<std::array<unsigned int, 2>, 18>
+ wedge_table_2{{{{0, 0}},
+ {{1, 0}},
+ {{2, 0}},
+ {{0, 1}},
+ {{1, 1}},
+ {{2, 1}},
+ {{3, 0}},
+ {{4, 0}},
+ {{5, 0}},
+ {{3, 1}},
+ {{4, 1}},
+ {{5, 1}},
+ {{0, 2}},
+ {{1, 2}},
+ {{2, 2}},
+ {{3, 2}},
+ {{4, 2}},
+ {{5, 2}}}};
+ } // namespace
+
+
+ template <int dim>
+ double
+ ScalarWedgePolynomial<dim>::compute_value(const unsigned int i,
+ const Point<dim> & p) const
+ {
+ AssertDimension(dim, 3);
+ AssertIndexRange(this->degree(), 3);
+
+ const auto pair = this->degree() == 1 ? wedge_table_1[i] : wedge_table_2[i];
+
+ const ScalarPolynomial<2> poly_tri(this->degree());
+ const Point<2> p_tri(p[0], p[1]);
+ const auto v_tri = poly_tri.compute_value(pair[0], p_tri);
+
+ const ScalarPolynomial<1> poly_line(this->degree());
+ const Point<1> p_line(p[2]);
+ const auto v_line = poly_line.compute_value(pair[1], p_line);
+
+ return v_tri * v_line;
+ }
+
+
+
+ template <int dim>
+ Tensor<1, dim>
+ ScalarWedgePolynomial<dim>::compute_grad(const unsigned int i,
+ const Point<dim> & p) const
+ {
+ AssertDimension(dim, 3);
+ AssertIndexRange(this->degree(), 3);
+
+ const auto pair = this->degree() == 1 ? wedge_table_1[i] : wedge_table_2[i];
+
+ const ScalarPolynomial<2> poly_tri(this->degree());
+ const Point<2> p_tri(p[0], p[1]);
+ const auto v_tri = poly_tri.compute_value(pair[0], p_tri);
+ const auto g_tri = poly_tri.compute_grad(pair[0], p_tri);
+
+ const ScalarPolynomial<1> poly_line(this->degree());
+ const Point<1> p_line(p[2]);
+ const auto v_line = poly_line.compute_value(pair[1], p_line);
+ const auto g_line = poly_line.compute_grad(pair[1], p_line);
+
+ Tensor<1, dim> grad;
+ grad[0] = g_tri[0] * v_line;
+ grad[1] = g_tri[1] * v_line;
+ grad[2] = v_tri * g_line[0];
+
+ return grad;
+ }
+
+
+
+ template <int dim>
+ Tensor<2, dim>
+ ScalarWedgePolynomial<dim>::compute_grad_grad(const unsigned int i,
+ const Point<dim> & p) const
+ {
+ (void)i;
+ (void)p;
+
+ Assert(false, ExcNotImplemented());
+ return Tensor<2, dim>();
+ }
+
+
+
+ template <int dim>
+ void
+ ScalarWedgePolynomial<dim>::evaluate(
+ const Point<dim> & unit_point,
+ std::vector<double> & values,
+ std::vector<Tensor<1, dim>> &grads,
+ std::vector<Tensor<2, dim>> &grad_grads,
+ std::vector<Tensor<3, dim>> &third_derivatives,
+ std::vector<Tensor<4, dim>> &fourth_derivatives) const
+ {
+ (void)grads;
+ (void)grad_grads;
+ (void)third_derivatives;
+ (void)fourth_derivatives;
+
+ if (values.size() == this->n())
+ for (unsigned int i = 0; i < this->n(); i++)
+ values[i] = compute_value(i, unit_point);
+
+ if (grads.size() == this->n())
+ for (unsigned int i = 0; i < this->n(); i++)
+ grads[i] = compute_grad(i, unit_point);
+ }
+
+
+
+ template <int dim>
+ Tensor<1, dim>
+ ScalarWedgePolynomial<dim>::compute_1st_derivative(const unsigned int i,
+ const Point<dim> & p) const
+ {
+ return compute_grad(i, p);
+ }
+
+
+
+ template <int dim>
+ Tensor<2, dim>
+ ScalarWedgePolynomial<dim>::compute_2nd_derivative(const unsigned int i,
+ const Point<dim> & p) const
+ {
+ (void)i;
+ (void)p;
+
+ Assert(false, ExcNotImplemented());
+
+ return {};
+ }
+
+
+
+ template <int dim>
+ Tensor<3, dim>
+ ScalarWedgePolynomial<dim>::compute_3rd_derivative(const unsigned int i,
+ const Point<dim> & p) const
+ {
+ (void)i;
+ (void)p;
+
+ Assert(false, ExcNotImplemented());
+
+ return {};
+ }
+
+
+
+ template <int dim>
+ Tensor<4, dim>
+ ScalarWedgePolynomial<dim>::compute_4th_derivative(const unsigned int i,
+ const Point<dim> & p) const
+ {
+ (void)i;
+ (void)p;
+
+ Assert(false, ExcNotImplemented());
+
+ return {};
+ }
+
+
+
+ template <int dim>
+ std::string
+ ScalarWedgePolynomial<dim>::name() const
+ {
+ return "ScalarWedgePolynomial";
+ }
+
+
+
+ template <int dim>
+ std::unique_ptr<ScalarPolynomialsBase<dim>>
+ ScalarWedgePolynomial<dim>::clone() const
+ {
+ return std::make_unique<ScalarWedgePolynomial<dim>>(*this);
+ }
+
+
+
template class ScalarPolynomial<1>;
template class ScalarPolynomial<2>;
template class ScalarPolynomial<3>;
+ template class ScalarWedgePolynomial<1>;
+ template class ScalarWedgePolynomial<2>;
+ template class ScalarWedgePolynomial<3>;
} // namespace Simplex
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2020 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE.md at
+// the top level directory of deal.II.
+//
+// ---------------------------------------------------------------------
+
+
+// Distribute Simplex::FE_Wedge on a DoFHandler.
+
+#include <deal.II/dofs/dof_handler.h>
+
+#include <deal.II/grid/grid_out.h>
+#include <deal.II/grid/tria.h>
+
+#include <deal.II/simplex/fe_lib.h>
+
+#include "../tests.h"
+
+using namespace dealii;
+
+void
+test_3()
+{
+ const int dim = 3;
+ const int spacedim = 3;
+
+ std::vector<Point<spacedim>> vertices;
+ std::vector<CellData<dim>> cells;
+
+#if true
+ vertices.emplace_back(0.0, 0.0, 0.0);
+ vertices.emplace_back(1.0, 0.0, 0.0);
+ vertices.emplace_back(0.0, 1.0, 0.0);
+ vertices.emplace_back(0.0, 0.0, 1.0);
+ vertices.emplace_back(1.0, 0.0, 1.0);
+ vertices.emplace_back(0.0, 1.0, 1.0);
+
+ {
+ CellData<dim> cell;
+ cell.vertices = {0, 1, 2, 3, 4, 5};
+ cells.push_back(cell);
+ }
+#else
+ vertices.emplace_back(0.0, 0.0, 0.0);
+ vertices.emplace_back(1.0, 0.0, 0.0);
+ vertices.emplace_back(0.0, 1.0, 0.0);
+ vertices.emplace_back(1.0, 1.0, 0.0);
+ vertices.emplace_back(0.0, 0.0, 1.0);
+ vertices.emplace_back(1.0, 0.0, 1.0);
+ vertices.emplace_back(0.0, 1.0, 1.0);
+ vertices.emplace_back(1.0, 1.0, 1.0);
+
+ {
+ CellData<dim> cell;
+ cell.vertices = {0, 1, 2, 4, 5, 6};
+ cells.push_back(cell);
+ }
+ {
+ CellData<dim> cell;
+ cell.vertices = {2, 1, 3, 6, 5, 7};
+ cells.push_back(cell);
+ }
+#endif
+
+ Triangulation<dim, spacedim> tria;
+ tria.create_triangulation(vertices, cells, SubCellData());
+
+ GridOut grid_out;
+ std::ofstream out("mesh.vtk");
+ grid_out.write_vtk(tria, out);
+
+ Simplex::FE_WedgeP<dim, spacedim> fe(2);
+
+ DoFHandler<dim, spacedim> dof_handler(tria);
+ dof_handler.distribute_dofs(fe);
+
+ deallog << dof_handler.n_dofs() << std::endl;
+
+ std::vector<types::global_dof_index> dof_indices;
+ for (const auto &cell : dof_handler.active_cell_iterators())
+ {
+ for (const auto face : cell->face_indices())
+ {
+#if false
+ dof_indices.resize(face <= 1 ? 3 : 4);
+#else
+ dof_indices.resize(fe.n_dofs_per_face(face));
+#endif
+ cell->face(face)->get_dof_indices(dof_indices);
+
+ for (const auto i : dof_indices)
+ deallog << i << " ";
+ deallog << std::endl;
+ }
+ }
+}
+
+int
+main()
+{
+ initlog();
+
+ test_3();
+}
--- /dev/null
+
+DEAL::18
+DEAL::1 0 2 6 8 7
+DEAL::3 4 5 9 10 11
+DEAL::0 1 3 4 12 13 6 9 15
+DEAL::1 2 4 5 13 14 7 10 16
+DEAL::2 0 5 3 14 12 11 8 17
\ No newline at end of file
projects / dealii.git / commitdiff